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‘Machine Proof of a Theorem on Cubic Residues 


By D. H. Lehmer, E. Lehmer, W. H. Mills, and J. L. Selfridge 
If p is a prime of the form 6m + 1, the numbers 
1°, 2°, ---, (p— 1)’ 
hen reduced modulo p consist of only (p — 1)/3 = 2m distinct numbers between 


land p — 1. These 2m numbers are known as the cubic residues of p. Thus, the 
pubic residues of 13 are 1, 5, 8, and 12, while those of 97, when arranged monotonely, 


1, 8, 12, 18, 19, 20, 22, 27, 28, 30, ---. 


Here we observe the triplet (18, 19, 20) of three consecutive numbers among the 
bic residues of 97, while no such phenomenon exists for p = 13. We will call any 
Sset of three consecutive positive integers a triplet. A prime p = 6m + | is called 
pptional if it does not have a triplet of cubic residues. Thus 13 is an exceptional 
ime, and 97 is not an exceptional prime.’ It has been known since 1928 [1] that 
“sufficiently large’ primes have a triplet of cubic residues. Thus there are only 


finite number of exceptional primes. By using machine methods we have proved 
much more, namely: 


THEOREM 1. 
(a) The only exceptional primes are 


2, 3, 7, 13, 19, 31, 37, 43, 61, 67, 79, 127, 283. 


(b) Every non-exceptional prime has a triplet of cubic residues that does not 
exceed 


(1) (23532, 23533, 23534). 


(c) There are infinitely many primes whose smallest triplet of cubic residues is 
(1). Hence, result (b) is the best possible. 

Remark. The referee comments that the proof of Theorem 1, described below, 
is “not a machine proof in the sense of the theorem-proving programs now being 
fleveloped.” This is true. The aim of most writers on this subject is to consider a 
very general program enabling a digital computer to prove a wide class of theorems 

t a very low level, beginning with the axioms, setting its own goals, and trying to 

hieve them without human intervention. This is, in a way, a simulation prob- 
em. Speculations about such programs involve (significantly) such notions as 
Mecidability. Meanwhile, no really new theorems seem to emerge. Perhaps too 
much is expected of a single program. 

In our work, instead of starting with axioms, we did not hesitate to use any 
Mevice or previously known result that might be useful. In particular, the authors 
aided and abetted the machine in its search for a theorem and its proof. Neverthe- 

Received April 3, 1962. 

1 Every prime p = 6m — 1 is non-exceptional since every number less than p is a cubic 


Tesidue of such a prime. 
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less, all three results (a), (b), and (c) are due to the machine. Even the verifica- 
tion of these results using the data supplied by the machine would be far too long 
and hazardous a calculation to do by hand. 

It is perhaps worth noting that (a), (b), and (c), though proved in a finite 
number of steps, are statements about infinite classes. For example (a) does not 
assert that the only exceptional primes less than one million are 7, 13, --- , 283. 
This would have been merely a new finite result, easily obtainable by the machine, 
but not a “genuine” theorem. 

We give in what follows an explan:.ion of how the computer was programmed 
to carry out the immense number of steps needed to prove this theorem. For dis- 
cussion of the general problem for runs of kth power residues, the reader may con- 
sult previous papers [3] and [4]. We note here that a corresponding theorem for 
pairs of consecutive cubic residues has been proved by M. Dunton [2], and that 
there is no such theorem for four consecutive cubic residues [3]. 

For a prime p = 6m + 1, the 2(p — 1)/3 = 4m non-residues fall into two 
classes such that the product of a cubic residue by a non-residue of one class is 
congruent modulo p to a non-residue of the same class. The product of two non- 
residues of the same class is congruent to a non-residue of the other class, while 
the product of two non-residues of different classes is congruent to a cubic residue 
of p. We call these two classes of non-residues Class 1 and Class 2 respectively. This 
definition becomes unambiguous as soon as one member is assigned to Class 1. 
Let Class 0 denote the class of cubic residues. Thus the numbers from 1 to p — 1 
are divided into three classes, each having (p — 1)/3 elements. We set R(s) = i, 
if s is congruent modulo p to a member of Class 7. Thus for every integer s not 
divisible by p, R(s) is defined and R(s) = 0, 1, or 2. Moreover it follows from 
the above discussion that 


(2) R(si82) = R(s:) + R(s2) (mod 3) 


for any s; and s2 not divisible by p. 
Next let S be a given finite set of distinct primes, say 


S={m,@,-°--, 4d, “<< >>-< &. 
The vector A = [a,, a2, --- , a], where each a; = 0, 1, or 2, will be called an 
S-vector. If p is a prime not in the set S, and R(qi) = a, R(m) = 
a@, +--+, R(q:) = a,, then A will be called an S-vector belonging to p. If all the 


primes q; in S are cubic residues of p, then the zero vector belongs to p. Except for 
this case, there are two S-vectors belonging to a given prime p, due to the choice 
in the definition of Class 1 and Class 2. 

There are 3‘ possible S-vectors. According to a theorem of Kummer (See [5], 
pp. 426-428), each of these 3‘ possible S-vectors belongs to an infinite number of 
primes. Thus the primes of the form 6m + 1, not in the set S, are divided into 
4(3° + 1) subsets, and each of these subsets contains an infinite number of primes. 

Now let n be an integer whose prime factors all belong to the set S, so that 


Nn=Qn'q@” ---q, (a; 2 0), 
and let 


a; = 3m; + 0b; (0 s b; 
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Thus 
(3) [b, , be, -++, bi 


is a sequence of ternary digits uniquely determined by n and S. We call (3) the 
decomposition vector of n. Suppose A belongs to the prime p. Then (2) yields 


R(n) = abi + abs + --- + ad, (mod 3). 
In particular n is a cubic residue of p if and only if 
(4) ab; + ab. + --- + ab, =0 (mod 3). 


More picturesquely we say that n is a cubic residue of p if and only if the decom- 
position vector (3) of n and an S-vector A belonging to p are mutually perpen- 
dicular or orthogonal modulo 3. 

If each member of the triplet (n — 1, n, n + 1) has its prime factors restricted 
to the set S, and if each of the three decomposition vectors of n — 1, n,n + 1 
respectively is orthogonal to the S-vector A modulo 3, then any prime p to which 
A belongs has this triplet of cubic residues. In this case we say that the triplet 
(n — 1, n, n + 1) disposes of the S-vector A. If we can dispose of each of the 
3° possible S-vectors this way using triplets not exceeding (VN — 1, N, N + 1), 
then it follows that every prime not in the set S has a triplet of cubic residues not 
exceeding (N — 1, N, N + 1). Finally, to show that there exist primes such that 
(N — 1, N, N + 1) is the smallest triplet of cubic residues, we must take S* to 
be the set of all primes less than N + 2 and find an S*-vector A* such 
that (V — 1, N, N + 1) is the smallest triplet that disposes of A*. If we can do 
all this with N = 23533, then we have proved Theorem 1. 

Before we made the machine runs we had no way of knowing that 23533 was 
the correct value of N, in fact we did not even know if there was a value of N with 
the required properties. Hence it was necessary to experiment with different values 
of N—in fact we made machine runs with seven distinct values of N. 

The success of our program depended very largely on the selection of a suitable 
set S. If po is an exceptional prime not in the set S, then no triplet will dispose of 
the S-vectors belonging to po . Hence if we are to dispose of all S-vectors, the set 
S must include all the exceptional primes. Unfortunately, while we know that the 
set of exceptional primes is finite, we have no way of finding them all in advance, 
or even of finding an upper bound for them. However a preliminary machine run 
was designed to test individual primes. With this program we tested all primes less 
than 11243. This gave us the exceptional primes 2, 3, 7, 13, 19, 31, 37, 43, 61, 67, 
79, 127, and 283. It seemed unlikely that there were any more, but we could not 
be sure of this until we ran the main program on the machine. 

The most important consideration in the choice of the set of primes S is the 
necessity of having a large number of suitable triplets available. We can estimate 
the number of triplets required by a crude probabilistic argument. If none of the 
three numbers, n -- 1, n, n + 1 is a cube and if all of their prime factors are in 
S, then the a priori probability that (n — 1, n, n + 1) disposes of a given S- 
vector A is 1/27. From this we estimate that the number of triplets required to 
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TABLE 1 
First Members of Triplets 





























t qt N; 
| 
Si, & 2 |1,2 
3 5 3 3, 4,8 
4 7 5 5, 6, 7, 14, 48 
5 1] 5 9, 10, 20, 54, 98 
6 13 9 11, 12, 13, 24, 25, 26, 63, 64, 350 
7 17 10 15, 16, 32, 33, 34, 49, 50, 119, 168, 440 
8 19 12 17, 18, 19, 38, 55, 75, 76, 152, 169, 208, 323, 2430 
9 23 11 21, 22, 23, 44, 68, 90, 160, 207, 322, 390, 2023 
10 29 17 | 27, 28, 56, 114, 115, 143, 174, 230, 288, 493, 550,782, 1274, 
2000, 3248, 9800, 13310 
11 31 16 29, 30, 31, 62, 91, 124, 153, 154, 340, 341, 494, 527, 713, 
1518, 1519, 13454 
12 37 24 35, 36, 37, 74, 110, 184, 185, 220, 259, 405, 406, 665, 702, 
960, 999, 1330, 1443, 1664, 1700, 2736, 3625, 5290, 7104, 
17575 
13 41 | 27 | 39, 40, 80, 123, 203, 245, 246, 285, 286, 287, 368, 492, 574, 
| 1023, 1024, 1188, 1517, 1680, 1681, 1885, 2294, 3772. 
| 4959, 5082, 29600, 32798, 212380 





For each t, g, denotes the tth prime and N, denotes the number of triplets in 
which the largest prime factor involved is q; . The smallest member of each of these 
N; triplets is listed. The fact that the table is complete is established in [6] 


dispose of all 3‘ vectors is approximately 


a 

log (27/26) 
However, the actual number of triplets available for small values of ¢ is much less 
than 29t. For example, if S is the set consisting of the first ¢ primes, then the num- 
ber of triplets < 442224 surpasses 29¢ only for ¢ = 25. For t = 13 the total number 
of such triplets, irrespective of size, is only 141. The first members of these 141 
triplets are given in Table 1. For t = 50, however, there are more than 1800 triplets 
less than 25000, which compares favorably with 29% = 1450. 

It is clear that before the proposed program can be attempted we must have 
the machine supply itself with a large lis. of triplets (n — 1, n,n + 1) whose prime 
factors are in S and where n S N. This preliminary program involves choosing a 
set S of primes and the limit N. For the first run one should take N rather large 
in order to be on the safe side. 

To decide quickly whether the prime factors of a number n S N belong to the 
set S there are three methods available. One may simply attempt to factor n 
using as trial divisors only the primes q; in S. If this fails to give a complete fac- 
torization, then n is not acceptable, since it contains a prime factor not included in 
set S. At least ¢ divisions are required for each nonacceptable n. A second method 
involves the construction of the product 


M = qi'qs? $019 at, 
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where §; is determined by the inequality 
ai < N < qft™. 


Then n is acceptable if and only if n divides M. The number M will ordinarily 
consist of several machine words. However, 


log M < tlog N 


so that approximately ¢ log N/log W division instructions need to be executed for 
each n, where W stands for the largest integer representable by a machine word. 
Since N < W, the number of divisions required is much less than in the first method. 
The third method consists of sifting out the multiples of primes not in the set S 
from the numbers from 1 to N. If a binary machine is used, a compact bit representa- 
tion technique is available in which one obtains a binary number of N bits whose 
nth bit is 1 or 0 according as n is acceptable or not. In this case it is particularly 
easy to look for three consecutive acceptable n’s that make a triplet. 

The problem of registering the list of triplets so as to make it most readily 
available to the main routine will be discussed later. We next consider the ordering 
of the list. 

If in the vector (3), bg + 0 and b, = 0 for all h > d, then we say that the 
vector is of dimension d. The dimension of the zero vector is defined to be zero. 
The dimension of a triplet (n — 1, n, n + 1) is defined to be the maximum of the 
dimensions of the three decomposition vectors of n — 1, n, n + 1 respectively. 
For example, if S contains the first 10 primes, then the triplet 


n — 1 = 9800 = 2°-5°-7? ~ (0, 0, 2, 2,0, --- , 0} 
n = 9801 = 3°-11° ~ (0, 1, 0, 0, 2, 0, --- , 0} 
n+ 1 = 9802 


2-13-29 ~ [1, 0, 0, 0, 0, 2, 0, 0, 0, 1, 0, --- , 0] 


has dimension 10. As a final step in our preparation of the list of triplets for the 
main routine, we sort the list by dimension and prepare a small table in which the 
machine can look up the address of the first triplet of each dimension d. 

In actual practice the number ¢ is large enough so that it is prohibitive to con- 
sider all 3‘ possible S-vectors separately. Therefore, it is necessary to dispose of 
many of them at once. We do this by means of the concept of case vector. Let 
d = t. A case vector of dimension d is a vector 


(5) C = [am ,a@,---, Gal, 


where a; = 0, 1, or 2. If we can find a triplet (n — 1, n, n + 1) of dimension d, 
whose three decomposition vectors are orthogonal to C modulo 3, then in one blow, 
we have disposed of all those 3°“ original S-vectors whose first d components agree 
with C. 

We can now describe the main routine and its methodical disposal of these case 
vectors. This is easily done with the flow chart in Figure 1 and a brief explana- 
tion. 





2 Figure 1 describes a program in which all 3‘ possible S-vectors are examined. However, 
by symmetry and other considerations, it is sufficient to examine the vectors from [0, 2] to [2]. 
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By “success” we mean, of course, that the machine has discovered a triplet of 
dimension d which disposes of the current case vector. As the proof proceeds the 
dimension d rises and falls, never rising more than a unit at a time, but often fall- 
ing more than a unit in an irregular way. In fact, if one imagines a ternary point 
to the left of a, , the case vector (5) becomes a rational variable, written to the 
base 3, that rises monotonely from zero to one with irregular speed, certain intervals 
being much more difficult than others. Since the dimension d rises by one unit each 
time there is no triplet available of dimension d which disposes of the current case 
vector, it follows that when the dimension becomes ¢ + 1, then the machine has 
examined the entire list of triplets without success. It then reports the S-vector 
as one that the data cannot handle. The processing of this output is discussed later. 

It is clear from the diagram that the machine is spending almost all its time 
examining the list of triplets. It is also clear that the basic operation here is that 
of finding the inner product of the case vector and a decomposition vector, so that 
every effort should be made to shorten this program loop. To this effect we exploit 
the fact that decomposition vectors are usually quite sparse, i.e., they consist mostly 
of zeros. Hence, it is advantageous to use a more condensed format in which only 
the nonzero components of the decomposition vectors are involved. In one such 
representation the coded word 


Fig. 1 


(6) U1, Vi 5 Us, V2 5 °** 5 Ue, mK 50 








let. of 
Is the 
1 fall- 
point 
o the 
rvals 
, each 
, case 
e has 
rector 
iater. 
time 
; that 
» that 
xploit 
10stly 
. only 
such 





MACHINE PROOF OF THEOREM ON CUBIC RESIDUES 413 


corresponds to the decomposition vector in which the (v; + 1)-st competent is 
u;, t* = 1(1)k, and all other components are zero. For example, the number 
n = 15678 = 2-3°-13-67 has the decomposition vector 


fi, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) 
which corresponds to the coded word 
1, 0; 2, 1; 1, 5; 1, 18; 0. 


Each u; is 1 or 2, and two binary bits are used to represent it. Suppose 2” < t < 2”. 
Then w bits are required for each v; , and the word (6) is composed of (w + 2)-bit 
subwords. 

With the decomposition vector written in the form (6), its inner product (4) 
with the case vector (5) is 


k 

» Ue; +41 - 

t=1 
Since we need only evaluate this modulo 3 we interpret u; = 1 as “add” and u; = 2 
as “subtract”. Thus, we compute the inner product without multiplications. The 
case vector (5) is stored in ¢ successive words in the memory, and the word (6) 
is unpacked by a series of left shifts in an obvious manner, so that each (w + 2)-bit 
subword constructs a command in which the operation is either addition or sub- 
traction, as determined by its first two bits, and whose address is determined by 
its remaining w bits. The last symbol, zero, in the word (6) is interpreted as “end 
of message,”’ and when that is reached the accumulator total is reduced to 0, 1, or 
2 by a few additions or subtractions of 3. A simple test for zero now suffices to tell 
whether n is a cubic residue of those primes p that correspond to the current case 
vector. 

We began with ¢ = 32, w = 5, so that w + 2 = 7, and with our 35-bit binary 
words we could use decomposition vectors with 5 or fewer nonzero entries. This 
was sufficient for every triplet that we had occasion to use. Later we increased ¢ 
to 50 and then to 55, so that we had to take w = 6, w + 2 = 8, and therefore 
we had to omit a few triplets which involve numbers with 5 distinct prime factors. 

In conclusion, we give a brief account of the successive runs made by the machine 
which culminated in Theorem 1. 

The preliminary search for exceptional primes had revealed only one of them 
greater than 127, namely 283. It was decided that for the first run the set S would 
consist of the first 31 primes, i.e., those primes £127, and 283. A high limit of 
N = 442224 was used, in the hope of providing the machine with an adequate 
supply of triplets. This actually gave us 1381 triplets. This run was successful in 
that all S-vectors were disposed of. At this point the machine had proved the (a) 
part of Theorem 1, so that we knew that our list of exceptional primes was com- 
plete. Furthermore, we knew that results of the type (b) and (c) were indeed pos- 
sible. 

This first run required 59 minutes on the IBM 704 during which time the 
machine considered about 101,000 case vectors. 

The limit was now lowered to N = 2” = 131,072 and a second successful run 
was made. The third run with N = 2" = 65536, however, resulted in the output of 
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eight vectors, which were easily disposed of by including the primes 137 and 139 
in the set S—an operation that did not require the use of the computer. In the 
fourth run with N = 2” = 32768 the machine reported 120 vectors, in two distinct 
regions, which it could not handle. We were able to dispose of these 120 vectors 
by increasing the set S. An attempt to make a run at 2° was abandoned because of 
the large amount of output due to the fact that our supply of triplets had fallen 
below 1000. It now became apparent that many more primes would have to be 
included in the set S. Accordingly ¢ was increased to 50 and a run was made with 
N = 3-2" = 24576. This resulted in an output of more than 512 vectors, all but 
4 of which started with 


(7) (0, 2, 0, 2, 0, 2, 2, 0, 2, 2, 2, 2, 2, 2, 2, 2, 0, -->]. 

By including the prime 271 all these vectors can be disposed of by the three triplets: 
8671 = 13-23-29, 8672 = 2°-271, 8673 = 3-7°-59 
9212 = ee 9213 = 3-37-83, 9214 = 2-17-271 
24388 = 2’-7-13-67, 24389 = 29°, 24390 = 2-3°-5-271. 


The remaining 4 vectors which started with 
[1, 2, 2, 2, 1,0, 1, 2, 1,0, 0, 1, ---] 
were all disposed of by the triplet 
20008 = 2°-41-61, 20009 = 11-17-107, 20010 = 2-3-5-23-29, 


which the machine did not possess because 20010 has five distinct prime factors. 

At this point we thought that Theorem 1 might hold with N = 24389 instead 
of 23533. A sixth machine run was made with N = 24389. There were a number of 
additional output vectors, but the extra vectors all started as in (7) and were dis- 
posed of with the prime 271 as before. We now had part (b) of the theorem with 
N = 24389. In an attempt to prove part (c) with this value of N, a “case test” 
program was written. In this program we let S* be the set of all primes less than 
24391, and we take a particular S* vector A*. The program then finds the smallest 
triplet that disposes of A*. We took for A* various extensions of the vector (7) 
in which nearly all the components corresponding to large primes g were 2. The 
largest triplet put out by these runs was 


23532 = 2°-3-37-53, 23533 = 101-233, 23534 = 2-7-41°. 


The vector A* which produced this result consists of 


0’s corresponding to g = 2, 5, 11, 19, 59, 79, 89, 113, 191, 211, 223, 229, 
269, 373, 577, 829, 839, 1613, 2393; 


1’s corresponding to g = 233, 313, 353, 919, 967, 2671 


(8) 


2’s corresponding to all other primes g < 23535. 


This and other outputs of the case test runs suggested that the primes q = 233, 
313, 331, 449, and 967 should be incorporated in the main run. The final successful 
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run was therefore made with N = 23533 and ¢t = 55. This run alone constitutes a 
proof of parts (a) and (b) of Theorem 1, while the case test of the vector (8) 
proves part (ce) of the theorem. 

The last main run was done on the IBM 7090 in about 40 minutes and required 
the examination by the machine of about 250,000 case vectors. 

The various machine runs were made at the Computer Centers of the Uni- 
versity of California at Berkeley, Los Angeles, and Livermore, and at the Uni- 
versity of Washington in Seattle, and we are grateful to the directors of these 
laboratories for the free use of their equipment. The machines used were the IBM 
701, 704, 709, and 7090. We are also grateful to John Brillhart, David Mapes, and 
Vance Vaughn for donating their time to this unsupported research. 
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Minimum Periods, Modulo p, of First-Order 
Bell Exponential Integers 


By Jack Levine and R. E. Dalton 


1. Introduction. The integers of the title, B(n), can be defined by the generating 
function, given by Bell [1, 2], 


(1.1) ct => Bin)=. 


n=0 


These numbers have been known for a long time and have a variety of interesting 
interpretations which include: 
(a) B(n) = the number of rhyming schemes in a stanza of n lines (attributed 
to Sylvester by Becker [3], 
(b) B(n) = the number of pattern sequences for words of n letters, as used in 
cryptology, Levine [4], 
(c) B(nm) = number of ways n unlike objects can be placed in 1, 2, 3, 
or 7 like boxes (allowing blank boxes), Whitworth [5, p. 88], 
(d) B(n) = number of ways a product of n (distinct) primes may be factored, 
Jordan [6, p. 179], Williams [7]. 
Epstein [8] extended the definition of B(n) to include all real and complex 
numbers n by means of the representation 


’ 


Ms 
|“ 


u 
So 
oo 


(1.2) B(n) = ar . 


He also gave several asymptotic formulas for B(n) in addition to the numerical 
values of B(n) for n = 1, --- , 20. This paper, as well as [2], contains numerous 
references dealing with these numbers. 

For computational purposes, various defining relations are known, for example, 


(1.3) Bin) = > a —— zG ) (r—k)’, 


r=1 k=0 


given by Bell [1], and Mendelsohn and Riordan [9]. This formula, (1.3), is equiva- 
lent to 


(1.4) B(n) = + S(n,r), 


where S(n, r) are Stirling numbers of the second kind, and which was obtained by 
Broggi [10] and Becker and Riordan [11]. Other references relative to (1.3) and 
(1.4) are found in Epstein [8]. 


(1.5) B(in+ 1) = (B+ 1)’, 


where on the right, B” is to be replaced by B(m) after expansion, was given by 
‘cogne [12]. (See also [1, 2, 11]). 
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The difference formula, 
(1.6) B(n) = A"B(1), 


Becker and Browne [3], was found to be the simplest for a digital computer, and was 
used in the computation of the B(n) given in the present paper. 
For a study of arithmetic properties of B(n), the congruence of Touchard [13], 


(1.7) B(n + p) = Bin) + B(n +1), = mod p, 
for p a prime, is basic. 


In addition, for our purposes, we mention the following congruence given by 
Hall [14], Touchard [13], and Williams {7}, 





(1.8) B(n + p”) = B(n+1) + mB(n), ~~ mod p. 
It is known that the (minimum) period of the sequence (reduced mod p) 
is a divisor of 
= 3. 
(1.10) ¥, =i 


and Williams [7] has shown this minimum period is precisely N, for p = 2, 3, 5. 

In this paper we extend these results to primes p > 5. The results obtained are 
stated in the theorem below. 

THEOREM. The minimum period, mod p, of the sequence B(QO), B(1), --- , B(n), 
of first-order Bell exponential integers is N, for p = 7, 11, 13, and 17. For the remain- 
ing primes p < 50, p = 19, 23, 29, 31, 37, 41, 42, 47, no known proper divisor, 
N, of N,, with N < 10° can be a period. 

In the course of the computations connected with this theorem the results of 
Cunningham [15] on factoring N, have been extended to include several new factors 
for certain p. These are exhibited in Table 3. 

In addition, the values of B(n), n < 74, have been computed, and are given in 
Table 1. This extends results of Gupta [16] for n < 50. Also, the values of B(n), 
mod p, (n S p, p < 50) are given in Table 2. Such values are needed in testing for 
periods. 


2. Computation of B(n). The symbolic binomial expansion (1.5), though useful 
in the computation of the first several B(n), becomes bulky and time-consuming as 
n increases, since each successive B(n) computed by this iterative scheme requires 
n — 2 multiplications and n additions involving larger numbers at each iteration. 
Formula (1.6), together with the initial values B(0) = 1, B(1) = 1, by contrast, 
requires but n — 1 additions for each new B(n). (See Becker and Browne [3)). 
Such a difference formula as (1.6) is ideally suited for a digital computer, since it 
substitutes fixed-point addition for multiplication in which accuracy to the unit’s 
digit must always be maintained. The only limitation which presented itself was 
the increasing size of the integers and differences involved. Using an octuple-precision 
addition subroutine, the numbers were generated on the difference table until a 
B(n) or a difference exceeded 80 digits, the capacity of a standard IBM card. This 
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TABLE 2 
B(n) mod p,0 Sn Sp, p S 50 











n 23 5 7 11 138 17 19 23 29 31 37 41 43 47 
0 ee ne es oe ele. ey gg 
1 Pee Se ee ee ee ee oe oe 
2 222 = = be Sa ae Se ae he oe ae ie 
3 4 2s ee oe ee Be Se ee Se ee ee 
4 014 215 15 165 16 6 16 1 1 1 
5 238 0 114 6 2 21 15 11 9 5 
6 05 8 1 13 19 O 17 18 39 31 15 
7 28 610 3 383 7 9 2 16 17 31 
8 4 6 9 17 O 22 17 33 4 12 4 
9 5 9 146 0 10 6 5 2 32 34 44 
10 221318 9 4 4 17 27 4 2% 
11 2915 4 1 @B ll 27 2 30 31 
12 11 11 5 20 138 15 O 27 27 O 
13 2 6 7 1 13 #1 3835 2 38 24 
14 BMwMB.7D@.45 1-58 
15 11 146 9 20 30 36 9 27 42 
16 14415 5 2 16 2 4 42 38 
17 2167 8 @ 19 19 12 
18 10 6 14 #3121 «WW 1 2 
19 2 9 20 21 28 16 20 44 
20 420 3 2 3 26 43 
21 16 15 25 32 22 27 5 
22 22 5 26 32 33 39 25 
23 225 19 6 3 36 29 
24 7 16 34 23 42 3 
25 24 2 0 3 27 20 
26 11 15 26 38 41 10 
27 21 16 19 13 42 O 
28 18 17 12 13 27 23 
29 2 12 21 35 1 30 
30 3 10 23 Il 22 
31 2 1 4 35 44 
32 35 2 21 18 
33 35 37 3 35 
34 -26 34 28 46 
35 i7 14 33 11 
36 6 35 32 37 
37 2 40 28 45 
38 31 3 25 
39 6 7 16 
40 5 12 17 
41 2 41 5 
42 17 13 
43 2 9 
44 23 
45 37 
46 19 
47 
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condition occurred during the computation of B(75). The program, which used the 
SOAP I assembly program, was used on an IBM 650 to compute the 75 numbers in 
73 minutes. A check was made with Gupta’s highest value, B(50), and the numbers 
were found to be identical. 


3. Factorization of N, , (p < 50). From a result of Fontene [17], it follows that 
all factors of N, are of the form 2kp + 1, when p is an odd prime. Using this in- 
formation, a program was developed for the Univac 1105 in the USE compiler 
language. This simply involved successive division of N, by divisors of the form 


P, = 2pk + 1, k = 1,2,3,--- 


until a zero remainder was reached. Since the routine was single-precision for 
the divisors, the P,’s were limited in magnitude to one accumulator length on the 
Univac 1105, or to values P, < 2”. 
Table 3 gives the N, and the factors thus obtained. 
The following is a summary of new prime factors and other information not 
contained in Cunningham [15, p. 72]. 
Case p = 17. N;- is completely factored into the three prime factors 10949, 
1749233, 2699538733. 
Case p = 19. No factors of Nig have been found, but Ny contains no factor 
<17,005,305 


TABLE 3 
N,’s and Prime Factors (Indicated by*) 








» £2 
Pp N; p—1 
5 N, = 781 = 11*-71* 
7 Ny, = 1 37257 = 29*-4733* 
11 Nn = 2 85311 67061 = 15797*-1806113* 
13. Nis = 2523 95922 16021 = 53*-264031*-1803647* 
17. Ni = 51702 51636 78960 47761 = 10949*-1749233* -2699538733* 
19 Nig = 1099 12203 09223 96438 40221 No known prime factors 
23 Ne = 94911 21818 11268 72883 43196 77753 = 461*-1289*- 
1597216194112486480522357 
29 Ne = 9 17030 76898 61468 33772 08150 52610 77188 02981 = 


59*- 16763* -84449* - 2428577 * - 14111459* -32037737880884399 

31 Nx = 56897 24710 24107 86528 70214 34301 97715 85348 24481 No 
known prime factors 

37 =Naz = 29 31981 93216 04953 92799 53613 49988 42485 03538 78009 
36166 51181 = 149*-1999*-7993*. (quotient > 40 digits) 

41 Ny = 33271 94076 58177 99967 83498 10240 83656 39964 72332 
54041 27485 81284 48841 = 83*. (quotient > 40 digits) 

43 Na = 4129 46984 92929 20838 07232 88782 88579 08531 14434 
61669 54570 31137 54094 99893 = 173*-6709*. (quotient > 40 
digits) 

47 Na = 84 30270 13796 61926 57970 97431 77268 05988 90944 54377 
04795 47313 54904 95405 42692 40497 = 1693*- (quotient > 
40 digits) 
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Case p = 23. No new prime factors of Nz have been found, but the third 
factor 1,597,216,194,112,486,480,522,357 contains no factor 
<59,929,399 

Case p = 29. Four new prime factors of Nx are 16763, 84449, 2428577, 
14111459. 


4. Determination of minimum periods, mod p. The knowledge of B(1), B(2), 
--+ , B(p), (or of any p consecutive B’s) will determine the complete set of B’s, 
mod p. Hence, if N be a factor of N, , to test for a period of the sequence {B(n)} 
mod 7, it is sufficient to calculate B(N + 1), B(N + 2), --- , B(N + p), mod p, 
and compare with B(1), B(2), --- , B(p), mod p. 

Furthermore, if N, can be expressed as a product of r factors, it is not necessary 
to test all possible combinations of factors for periods, but merely the combinations 
of r — 1 factors. A positive result would indicate what further testings are necessary. 

In case the complete factorization of N, into prime factors is unknown it may 
not be possible to find the minimum period. 

The actual testing of the various factors for the period property was accom- 
plished on an IBM 650. The program requires N, the factor to be tested; p, the 
particular prime; and B(1), B(2), --- , B(p), mod p. These B’s were obtained from 
a modification of the program used to calculate Table 1 and are given in Table 2. 
The program used could test any factor less than 10“. It would, of course, be im- 
practical to calculate every B through B(N + p), so a process of proceeding in 
jumps of powers of p by means of (1.8) is used. 

The factor N being tested is first expressed to the base p, 


(4.1) N = a,p" + aap” + --- + ap+a. 


The various steps are then (all calculations mod p): 
(1) Calculate B(p + 1) by (1.7). 


(2) Calculate B(a,p” + x), x = 1, 2, --- , p + 1, by the iterations 
(4.2) B(tp"+y) = B(é-— 1p" + y+ 1) + nmB((E— 1p" + y), 
(4.3) B(tp" + p +1) = B(tp” + 1) + Bitp” + 2), 
where ¢t = 1, 2, --- , an; y = 1, 2, --- , p. Equation (4.2) follows from (1.8), and 
(4.3) from (1.7). 
(3) Calculate B(anp” + a,sp” + 2), x = 1, 2, ---,p +1, by 
(4.4) B(up"* +z) = B((u—1)p"* +241) + (n— 1)B((u—1)p""* +2), 
(4.5) B(up™* + p + 1) = B(up"™* + 1) + B(up™* + 2), 
where u = 1, 2, --- , Gnt;2 = Gap” + 1, --+, Gap” + p. 
This procedure is continued until we reach 
(4.6) B(M + 1), B(M + 2), ---, B(M + p), 
where 


M = asp” + anap” + --- + ap. 
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Since one member of (4.6) is B(N), we start from that point and calculate 
which are then compared with 


B(1), B(2), ahr , B(p), 


for the period property. The results of these calculations have been given in the 
theorem of Section 1. 
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Concerning the Numbers 27° + 1, p Prime 
By John Brillhart 


1. Introduction. In a recent investigation [7] the problem of factoring numbers 
of the form 2°” + 1, p a prime, was encountered. Since 2”? + 1 = (2? — 2*?*” +1) 
(2? + 2*?*» + 1) for odd p, the problem consists of factoring the two trinomials on 
the right. In this paper the results of a search for factors of these trinomials are 
given, as well as a determination of the nature of certain of these numbers for 
which no factor was found. 


2. Elementary factors. Let N, = (2? — 2**? + 1) (27 + 2H? 4.1) = 4, 
-B, , p an odd prime. 

A. From the fact that 5 | N,, it easily follows that 5 | A,iffp = +1 (mod 8) 
and 5 | B, iff p = +3 (mod 8). On the other hand, 5° + N, unless p = 5; for, 
since 2 is a primitive root of 25, 2 belongs to the exponent (25) = 20. But 2”? = —1 
(mod 25), or 2° = (mod 25). Therefore, 20 | 4p, or p = 5. Thus, if p = 5, 
5° | 2° + 1 = 1025, while if p ¥ 5, 5° + N,. 

B. If q is a prime “5 and q | N,, then 2° = 1 (mod q). But then 2 belongs 
to the exponent 4p (mod q). Thus by Fermat’s Theorem, 4p | g — 1; that is, every 
prime divisor #5 of A, or B,is =1 (mod 4p). 

C. Suppose p is odd and q = 4p + 1 isa prime. Then 2° = 2° =1 (mod q). 


It follows from Euler’s Criterion that 2” = (2) (mod q). But since p is odd, 
q@=5 (mod 8). Therefore, 2°” = —1 (mod gq), or q| 2°” + 1. Unfortunately, 


however, it has not been possible to discover the conditions that determine which 
of A, and B, q will divide. 


3. The Search. 

A. Extent. The search for prime factors g ~ 5 of A, and B, , which was conducted 
on the IBM 701 at the University of California, Berkeley, was made over the follow- 
ing intervals: 


1<q<~WJ/Bws for By 
1<q<3-2” for An 
Yeect for 71 < p < 179and p = 241 
i<eet for 179 < p < 1200, p ¥ 241. 


No N, for p < 71, p + 59, were considered, since these numbers have been 
completely factored. Nos was examined along with N;; to the bound 2”, these num- 
bers being of particular interest (See [7]). 

B. Results. (i) The program produced a vast number of new factors, as well as 
several corrections to the literature (See [4]). The new factors of N,, p < 250, are 
indicated in the accompanying table by * to distinguish them from factors pre- 
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viously known [2]. For 250 < p < 1200 all factors > 300,000 are new, and are 
therefore not indicated by *. A dot following the final factor means that the nature 
of the complementary factor is unknown. 

(ii) A complete factorization was accomplished for By, As; , and Aj, the 
primality of the complementary factor in each case being assured by the non-ex- 
istence of a factor below its square root. The factorization of By is of particular 
interest, since this number appears in [2] and [3] as a prime. 

The author would like to thank Mr. K. R. Isemonger for providing the com- 
plete factorization of By;, as well as the much sought after factorization for An , 
which, previous to his attack on the number, had only been known to factor into 
the product of two primes. 

(iii) A program was written to test the divisibility and multiplicity of all known 
factors, with the result that all factors were found to be correct, but none was found 
to be multiple. 

C. The Program. The structure of the search program was similar to that de- 
scribed in [1]. In particular, for each p a table of differences was computed from the 
first 1155 = 3-5-7-11 terms of the sequence 4pk + 1, k = 1, 2, ---, that re- 
mained after the multiples of 3, 5, 7, and 11 had been sieved out. This table was 
used repeatedly by the program to produce a sequence of trial divisors, among 
which the factors, if any, were to be found. The remainders of A, and B, for each 
trial divisor were calculated by residue methods, both remainders being calculated 
at the same time because of the similarity in form of A, and B, . The occurrence 
of a 0 remainder in this calculation signalled the discovery of a factor of one of 
the two numbers, but not both, since obviously they are relatively prime. To ex- 
amine each NV, required from 5 to 15 minutes, the N, for the larger p’s requiring a 
shorter time. 


4. Primality Testing. 

A. At the conclusion of the search for factors, the primality of several numbers 
of immediate interest, namely, Az; and Ao , was still in doubt, because no factor 
had been found. It was then noted by Professor D. H. Lehmer that the primality 
of numbers of the form under consideration could be decided by Proth’s Theorem 


[5]: “If M = k-2” + 1, where 0 < k < 2”, and (*) = — I, then M is prime iff 
a’““- = —1 (mod M).” In the present case A,, B, = M = (2*”"” + 1)-240* 


+ 1, withO < k = 2?" +1 < 2*°* for p an odd prime, the value of a being 
easily obtained from the reciprocity law for the Jacobi symbol. 

A program was accordingly written by Professor Lehmer for the IBM 701 to 
calculate the required residues. The modulus used for each test was N, rather than 
the A, or B, in question, so that the reduction of the successive powers could be 
accomplished by multi-precision subtraction instead of division by a multi-precision 
divisor. The remainder thus produced was further reduced mod A, or B, by a sub- 
tractive routine written by the author. The final residues in binary from both rou- 
tines have been preserved on IBM cards for later checking purposes. 

B. It is believed that the two testing programs were accurate, since the antici- 
pated results were obtained in every trial case save one. In this case, By , a discrep- 
ancy existed between the literature, which stated the number was prime, and the 
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TABLE OF Factors 








p Qr — Qhe+tD 4 4 Qp 4+ Qurt) 4 
3 5 13 
5 5? 41 
7 113 5-29 
11 5-397 2113 
13 5-1613 53-157 
17 137 -953 5-26317 
19 5-229 -457 525313 
23 277 - 30269 5- 1013-1657 
29 5- 107367629 53690368 1 
31 5581 - 384773 5-8681-49477 
37 5-149- 184481113 593 -231769777 
41 181549 - 12112549 5- 10169 -43249589 
43 5- 1759217765581 173 - 101653 - 500177 
47 140737471578113 5-3761- 7484047069 
53 5-1801439824104653 15358129 - 586477649 
59 5- 1181-3541 -157649- 5521693* - 104399276341 * 
174877 
61 5-733 - 1709 - 368140581013 3456749 - 667055378149 
67 5-269 -42875177 - 15152453 - 9739278030221 
2559066073 
71 4999465853 - 472287102421 5-569 - 148587949 - 5585522857 
73 prime 5-293 -9929 -649301712182209 
79 prime 5-317- 
83 5- 13063537 *- 997 - 
148067 197374074653* 
89 1069- 5: 
97 389 - 4657 - 5- 3881-5821 - 3555339061 - 
394563864677 
161 5- 809- 
103 41201 - 520379897 *- 5-17325013*- 
4730001577 11296729* 
107 5-857- 843589 - 
109 5- 5669 - 666184021 *- 
113 prime 5- 58309 -2362153*- 
127 509 - 26417 - 140385293*- 5- 18797 -72118729*- 
131 5-642811237*- 269665073* - 
137 189061 - 5- 
139 5-1408349*- 557 - 
149 5- 1789- 
151 prime 5: 
157 . prime 
163 5-653 -9781 -7807049*- prime 
167 prime 5-75005713*- 
173 . c 
179 5-31815461*- c 
181 5-9413- c 
191 25212001 *- 5-3821- 
193 773: 5-3089 - 148997 - 
197 5-4729- 52009 - 
199 797 - 5: 
211 5-95110361*- c 
223 95768689* - 5-11597 -6530333*- 
227 5- 54449 -83132849*- 




















CONCERNING THE NUMBERS 2” + 1, p PRIME 


TABLE OF Facrors—Continued 


427 




















p 27 — Quet+) 4 1 2° 4+ Quet) 4 1 

229 5- 2749 -5523481*- c 

233 30757 - 5-3108221*- 

239 prime 5- 

241 prime 5- 2640397 * - 15594629* - 

251 5-1912621- 5021- 

257 c 5- 28564009 - 

263 c 5-119929-731141- 

269 5-2153 - 3229 - 5381 - 8609 - 
4273873 - 

271 10474693 - 5-97561- 

277 5-1109- 23268 1 - 98002601 - 

281 91568909 - 5-3373 -3827221- 

283 5- prime 

293 5- 22396921 - 5861 - 12893 - 60488093 - 

307 5-93329 - 1021697 - 1229 -7369 - 254197 - 201846361 - 

311 6221 -21149- 5- 

313 42569 - 68 1089 - 6386453 - 5- 

317 5: c 

331 5-589181- c 

337 683437 - 30499849 - 5- 5393 - 32353 - 2549069 - 

347 55575597 - 60988721 - 2777 - 

349 5-8377 - 763613 - c 

353 prime 5- 

359 585889 - 5199757 - 5- 

367 prime | 5- 

373 5-1493- | ¢ 

379 5-4549 - 10219357 - | prime 

383 13789 - 111650629 - 5-4597 

389 5-17117-51349 -2852149- c 

397 5-11117- 14293 - 25409 - 6312301 - 

401 c 5-3209- 

409 1637 -9817- 5-4909 - 1531297 - 1856861 - 

419 5- 63689 - 356989 - 53633 - 186037 - 

421 5-31142213- c 

431 91373 -3754873- 5: 

433 1733 -5197- 5-31177 -239017 - 

439 695377 - 5: 

443 ¢ 5- 

449 3615349 - 111190361 - 5-3593 - 165233 - 

457 rime 5-71293- 

461 5- 14753 - 7278269 - 226813 - 21102737 - 

463 c 5-46475941 - 

467 5- 13453337 - 252181 - 1372981 - 

479 6380281 - 39557737 - 5- 70309537 - 
79190197 - 

487 1949- 5-7793 -890237 - 

491 5- 3929 -34631213- 

499 5-43913 - 1179637 - 1997 - 

503 6037 - 10061 - 5- 

509 5- 103837 - 4073 - 13350053 - 

521 c 5-16673- 

523 5-8369 - 351457 - : 
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p Qr — QhetD 4 4 2? + Quet) 4 4 
541 5- 1281089 - 10393693 - 262302769 - 
547 5-67887077 - c 
557 5- c 
563 5- 51797 - 133489553 - 
569 37690561 - 5-47797 - 170701 -257189- 
571 5 - 2384497 - 5536417 - c 
94600997 - 
577 2309 - 92936237 - 5- 
587 5-35221- 13658317 - 
593 c 5- 
599 306689 - 9385133 - 5-4793 -86257 - 
601 7213: 5-79333 -685141 - 
607 c 5- 
613 5- 17458241 - 
617 c 5-86381- 
619 5- 114519953 - 2477 - 103993 - 284741 - 
631 ce 5-328121-651193- 
641 ce 5 - 62248793 - 
643 5- c 
647 144563093 - 5-854041 -9679121- 
653 5- c 
659 5§-5273 1534153 - 
661 5- c 
673 2693 - 26921 - 419953 - 5- 
4118761- 
677 5-5417 c 
683 5- ce 
691 5- 11057 - 
701 5- c 
709 5- 2837 - 
719 c 5-8629- 
727 2909 5- 
733 5- 627449 - 
739 5§- 523213 - 170756297 - 2957 -6139613- 
743 260683037 - 5- 
751 c 5-9013- 
757 5- c 
761 82189 - 529657 - 1567661 - 5-9133- 
769 c 5- 
773 5-9277 -961613 -8979169- 
28764877 - 
787 5- 47221 -406093 - 14121929. 
797 5- 
809 5- 6473 - 25889 - 1948073 - 
811 5- 5336381 - 
821 5- 
823 19753 - 17678041 - 5- 
827 5- 36389 - 148861 - 2312293 - 
829 5- 
839 5564249 - 5- 
853 5-3413- 


857 
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p 2? — 2h) 4+ 1 2? + 2QhetD 4 1 
859 5-82488053 - 41233 - 18970157 - 
863 62137- 5- 

877 5- 136813 - 178909 - 

881 292493 - 5- 

883 5- 3533 - 10597 - 

887 5- 

907 5- 

911 109321 - 5-29153- 

919 15174529 - 5- 3677 - 169097 - 
929 11149-319577 - 5-7433 -85469 -858397 - 
937 5-802073 - 

941 5- 3383837 - 

947 5- 189401 - 6522937 - 

953 5- 

967 328781 - 12056557 - 5- 47054221 - 
971 5: 19421 - 

977 5: 

983 5: 

991 47569 - 5-27749- 

997 5 - 3989 - 23929 - 1316041 - 

1009 12109- 5-242161- 

1013 5- 33449261 - 

1019 5-61141 - 207877 - 

1021 5- 88557457 - 

1031 181457 - 5 -32993 - 

1033 5-4133 -785089 - 
1039 4157 -47577889 - 5: 

1049 4640777 - 5- 

1051 5 -92489 - 2030533 - 1513441 - 77933753 - 
1061 5 -49459577 - 

1063 4253 - 119057 - 2351357 - 5: 

1069 5+ 25657 - 

1087 5-4349 - 182617 - 
1091 5-13093- 

1093 5- 13155349 - 4373 - 

1097 5- 114089 - 79321877 - 
1103 132361 - 5-525029 - 

1109 5-13309- 115337 - 

1117 5-67021 - 40213 -715148089 - 
1123 5-40429 - 4493 - 597437 - 
1129 5-4517- 

1151 5+ 36833 - 

1153 152197 -67796401 - 5: 

1163 5-37217 -37453253 - 

1171 5- 13152673 - 

1181 5: 1369961 - 9178733 - 

1187 5-9497 - 151937 - 

1193 3: 
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test routine, which stated the opposite. The number was immediately run on the 
factoring program, and rauch to the satisfaction of all concerned, a factor was found, 
and the test routine was exonerated. 

A further verification of a kind has come from Mr. Isemonger, who, acting on 
the test results that A; and By; were composite, succeeded in finding the factoriza- 
tions mentioned above. 

C. All A, and B,, 71 < p S 757, for which no elementary or other factor was 
known, were tested for primality. In all, 50 numbers were tested, with the result 
that 14 of them were found to be prime. These are listed as prime in the accompany- 
ing table, while the remaining 36 composite numbers are indicated as such by a 
““e” in the proper positions of the table. 

Each number with 71 < p < 457 was tested twice with complete agreement in 
the results. No number for p > 457 was tested twice, for testing a single number in 
this range required approximately 30 minutes. 
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Runge-Kutta Methods with Minimum 
Error Bounds 


By Anthony Ralston 


1. Introduction. Numerical methods for the solution of ordinary differential 
equations may be put in two categories—numerical integration (e.g., predictor- 
corrector) methods and Runge-Kutta methods. The advantages of the latter are 
that they are self-starting and easy to program for digital computers but neither 
of these reasons is very compelling when library subroutines can be written to 
handle systems of ordinary differential equations. Thus, the greater accuracy and 
the error-estimating ability of predictor-corrector methods make them desirable 
for systems of any complexity. However, when predictor-corrector methods are 
used, Runge-Kutta methods still find application in starting the computation and 
in changing the interval of integration. 

If, then, Runge-Kutta methods are considered in the context of using them for 
starting and for changing the interval, matters such as stability [2], [3] and minimiza- 
tion of roundoff errors [4] are not significant. Also, simplifying the coefficients so 
that the computation will be speeded up is not important and, on modern computers, 
minimization of storage [4] is seldom important. In fact, the only criterion of sig- 
nificance in judging Runge-Kutta methods in this context is minimization of 
truncation error. It is the purpose of this paper to derive Runge-Kutta methods of 
second, third and fourth order which have minimum truncation error bounds of a 
specified type. We will consider only the case of integrating a single first-order 
differential equation because this is the only tractable case analytically. But it 
seems reasonable to assume that methods which are best in a truncation error 
sense for one equation will be at least nearly best for systems of equations. 


2. The General Equations. For the solution of the equation 


(2.1) y’ =f(z,y) (20) = Yo 
at a sequence of points x; , 12, ---- the general Runge-Kutta method of order m is 
(2.2) Yaoi — Ya = k = wk; 


t=1 
where y, = y(2,), the w,’s are constants and 
i-1 
(2.3) k; = haf 2 + ajh, > Un + > But; 
j= 
with ha = 2n41 — Zn and a = 0. By choosing the a,’s and §;;’s properly the ex- 


pansion of the right hand side of equation (2.2) about (2, , y,) in powers of h, can 
be made identical with the Taylor series expansion of k about z, through the term 
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in h,”. One set of constraints on the a;’s and 6;;’s that always results is 


i—1 
(2.4) a= D> Bi; 

The error to be added to the right hand side of equation (2.2) to make it an 
exact relationship consists of a term Eh7**, where E depends on f(z, y), plus terms 
of higher degree in h, . For the classical fourth-order method of Kutta [6] a bound 
on E has been found by Lotkin [7] who improved on a bound of Bieberbach [1]. We 
will derive bounds on E for Runge-Kutta methods of second, third and fourth order 
in the same form as Lotkin: 


(2.5) | E| <cML” 

where c is a constant and, in a region R about (2, , yn)* 

(2.6) | f(z,y)| <M and te L***/M* 
: ; axtay? 


where M and L are constants and i + j < m. 


3. Second-Order Systems. The coefficients to be determined are w,, wz, 
oy and 62 . Matching powers of h, through h,” imposes three constraints (see, for 
example, (5]) which leads to the following one-parameter family: 


(3.1) w; = 1 — 1/(2a2) we = 1/(2a2) 

with Bx, = a» following from equation (2.4). The coefficient c; of h,’ is given by 
(3.2) ce = [() — (as'we/2)|D*f + (8)f, Df 

where 

(3.3) D = 8/dx + fr(9/dy) fn = f(Xn, Yn). 

Using (3.1) and (3.3) in (3.2) and the notation of (2.6), a bound on £ is given by 
(3.4) | EB | <[4| (8) — (a2/4) | + 3)ML* 


Clearly the minimum bound will be achieved if we set a2 = } which leads to the 
well known formula [5] 


(3.5) Yast — Yn = Bhaf(tn, Yn) + (¥)Rnf(tn + (F)hn , Yn + (F)hnfe) 
for which the right-hand side of (3.4) becomes (4) ML’. 
4. Third-Order Systems. In this case the coefficients form a two-parameter 
family given by 
w, = 1 + [2 — 3(a2 + as)]/6oras 
w2 = (3a; — 2)/[6a2(a3 — a2)) 
Ws = (2 — 3a2)/[603(a3 — a2)] 
Bsz = [a3(a3 — a2)]/[o2(2 — 3a2)] 


* This region must, of course, include all values of z and y in equations (2.2) and (2.3). 


(4.1) 
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when a: ~ 0, a3 ~ 0, and a2 ¥ a3. When a; = 0 there is a one-parameter family 
given by 


oe = w= }-—w 
(4.2) es i 
Bae 1/(4w3) Ww: = §. 
When a. = a; = 3 there is a one-parameter family given by 


(4.3) ws = we=i—ws Ba = 1/(40) 


No solutions other than (4.1), (4.2), and (4.3) are possible. 
The coefficient c; of ha’ is given by 


¢s = [(1/4!) — (1/3!) (exw. + a5'w;)] D*f 
(4.4) + [(1/4!) — (1/2!)ax’Baws)f, D*f 
+ [(3/4!) — arcsBx0s) Df Df, + (1/4)f,° Df. 


Substituting (3.3) and (4.1) into (4.4), combining terms, taking absolute values, 
and using (2.6) leads to the following bound on E: 


(4.5) |E| <[8|a:| +] a2| + | 2a2+ a3| + | a2 + a3| +2|as| +2] as |JML’ 
where 

a; = (xe) — (xe) [2(a2 + a) — Sara] 
(4.6) a2 = (gy) — 2/12 

a3= (§)—ax/6 aa = XX. 


A simple analysis shows that the coefficient of ML’ in (4.5) will be minimized if 
a2 = 4 and a; = } in which case we get 


(4.7) |E| < (})ML’. 
If (4.2) instead of (4.1) is used in (4.4) the bound on £ is 
(4.8) |E| < (#)ML’ 


independent of the value of w; . If (4.3) is used the bound, again independent of 


W3 , 18 
(4.9) |E| < (4)ML’. 


Thus, for third-order Runge-Kutta methods the minimum error bound of the type 
we are considering is given by (4.7). In this case equations (2.2) and (2.3) become 


(4.10) Ynsi — Yn = ($)Ki + (4)he + ($)ks 
where 

ky = Wnf(Zn , Yn) 
(4.11) ke = Raf(Zn + dhn , Yn + $h1) 


ks 


hif (Tn + (3)hn » Yn + ( $)ke). 











434 ANTHONY RALSTON 
5. Fourth-Order Systems. Again, the coefficients form a two-parameter family. 

They are given by 

w, = $+ [1 — 2(a2 + az)]/12a20;3 

We = (2az — 1)/[12a2(a3 — a2)(1 — ae)] 

ws = (1 — 2az)/[12a3(a3 — a2)(1 — as)] 
(5.1) ws = 3 + [2(a2 + a3) — 3]/[12(1 — a2)(1 — as)] 

a= 1 Bs = [as(as — az)]/[2a2(1 — 2ez)] 


(1 — a)[ar + as — 1 — (2a; — 1)’) 
2a2( as ~ a2) [Baz as —_ 4( a2 + a3) + 3] 


(1 — 2ae)(1 — a2)(1 — az) 
a3(a3 — a2)[6a2a3 — 4(a2 + a3) + 3} 


when az ~ 0, a3 ~ 0, a2 ¥ 1, as ¥ 1, and a, ¥ a;. The other possible solutions are 


Be = 








Baz = 


a2 = a, = 3 a = 1 
(5.2) Ws >= 4 WwW. = (2) — Ws w= 4 
Bs = 1/(6ws) Be = 1 — 3; Bag = 30s; 
and 
eg = 1 => 5 o- 1 
(5.3) w= % We = % — WH ws = 2 
Bse = § Be = — 1/(12w,s) Bas = 1/ (31) 
and 
a2 = } a; = 0 a = 1 
(5.4) w=(3)—-w w=} w=} 
Bs: = 1/(12ws) Bae = 3 Bas = 6; . 


For a. = 0 and a; = 1 no solutions are possible. Proceeding as before, a tedious 
computation leads to a bound on the error of 


| E| < [16| bi: | + 4] b2| + | be + 3b3| + | 2b, + 3b3| + | be + bs | 
(5.5) + |bs| + 8|bs| + | bs| + | 2bs + by[ + | bs + be + bz| 
+ | be | + | 2be + br| + | br| + 2] bs| JML* 


where 
bi = (ria) —- (ae) (cow, +> as Ws + ws) 
be = (ay) — Powers Baws + (cxBe2 + a18ss) i) 
bs = (rho) — (8) [o2'Bsews + (a2'Ber + o5'Bes) ws] 
(5.6) bs = (a5) — Flax cxBss + (a2'Ber + a3 Bis) wal 
bs = (rho) — For BsoBasto 
be = (as) — Ffax'Bso'ws + (cBee + as8es)’wr] 
br = (ria) — a2(1 + as) Bs2Bagws bs = riv- 
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Finding those values of a2 and a; which, after substitution of (5.1) into (5.6), mini- 
mize the coefficient of ML‘ in (5.5) is extremely tedious but not impossible. The 
result is ag = .4 and a; = (4) — (#5)(5)"” and in this case (5.5) becomes 


(5.7) |E| < 5.46 X 10°ML*. 

If (5.2) is used instead of (5.1), the bound is minimized if w; = $, in which case 
(5.8) |E| < (a8s)ML‘ = 7.22 x 10°M LL. 

Using (5.3) the bound is minimized if w, = 4%, in which case 

(5.9) |E| < (ges) ML‘ = 19.72 X 10°M LL. 

Using (5.4) the bound is minimized when w; = — + in which case 

(5.10) |E| < (44¢)ML' = 17.64 x 10° ML. 


Thus, the best bound is given by (5.7) and the complete set of equations in this 
case (correct to eight decimal places) is 


(5.11) yaar — Yn = -17476028k, — .55148066K, + 1.20553560k, + .17118478k, 


where 


ky = haf(Xn , Yn) 

ke = Iaf(tn + -4ltn 5 Yn + Aki) 

haf (Zn + .45573725hn , Yn + -29697761k; + .15875964he) 

ks = Iaf(tn + hn, Yn + .21810040k, — 3.05096516k2 + 3.83286476ks). 


(5.12) 


v3 


The coefficients for the classical fourth-order method of Kutta are given by (5.2) 
with w; = }. For this case Lotkin [7] found the bound 


(5.13) |E| < (f3s)ML‘ = 10.14 X 10°ML' 


which has a coefficient almost twice as great as that given by (5.7). This classical 
method may be considered a special case of the class of methods in which a, = 1 — 
a2 . The method of this class for a2 # 4 which has the minimum error bound also 
has reasonably simple coefficients and, therefore, is of some interest since it provides 
an improvement of the classical method even for hand computation. This method 
also has ag = .4 and the complete set of equations is 


(5.14) Ynsi — Yn = (ee) (11k, + 25ke + 25ks + 11k,) 
where 

ky = haf(Xn , Yn) 

ke = Nnf(tn + (¥)hn , Yn + ($)Ki) 


(5.15) 


> 
i 
| 


= Anf(tn + (¥)ha, Yn — (Po)hi + ($e) 

ky = Iaf(tn + hn, Yn + (ae) (19k, — 15k: + 40k) ). 
The error bound is given by 
(5.16) |E| < (aes) ML‘ = 7.70 X 10°M LL". 
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For completeness we list here the error bounds for two other fourth-order Runge- 
Kutta procedures which have appeared prominently in the literature: 
(1) a2 = }, a3 = 3 due to Kutta [6]. 


(5.17) |E| < (is) ML‘ = 9.91 X 10°M 
(2) a2 = 4, a3 = 3, ws = 1+ 1/(2)"” due to Gill [4]. 
(5.18) |E| < [(s8s) — (2s) 1ML' = 8.83 x 10°ML*. 


6. Numerical Examples and Conclusions. It is standard procedure in papers 
of this kind to finish with some numerical examples that illustrate how favorably 
the methods derived compare with other methods. Milne [8] remarks that especially 
in papers dealing with Runge-Kutta methods examples tend to be chosen which 
favorably illustrate the derived method. In fact, it is difficult to choose meaningful 
examples to illustrate Runge-Kutta methods and the reason for this is clear; the 
complicated nature of the error term makes it difficult to choose functions f(z, y) 


TABLE 1 
Errors in Integration of y’ = 1 — y*, y(0) = 0 (Solution: y = tanh z) 





| Magnitude of Error 
Step Size (h) | Number of Steps | Method | after Number of Steps 


| in Column 2 (X 10°) 
| 





12 
34 
65 


75 
138 
152 


1190 
2061 
2492 





low> OP a> 


| 
| 





TABLE 2 


Errors in Integration of y' = [e(y* + zy*® + 1)]/[8y°(zer — 6)], 
y(0) = 1 (Solution: y = [(@ + 5)/(6 — ze*)}"* 








| | Magnitude of Error 
Step Size (h) | Number of Steps | Method | after Number of Steps 


in Column 2 (X 10°) 








A 5 

a ) B 10 
| C 8 

A 26 

1 10 B 6 
C 14 
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which really serve as a test and at the same time are such that the differential «qua- 
tion (2.1) can be solved analytically. 

Here, for the sake of illustration and to indicate that the derivations in this 
paper were performed correctly, we present in Tables 1 and 2 the results of two 
simple computations comparing three Runge-Kutta methods:* 

Method A—Equations (5.11) and (5.12) 

Method B—Equations (5.14) and (5.15) 

Method C—The classical method with coefficients given by (5.2) with w; = }. 
In the example in Table 1 method A compares quite favorably with the other two 
while the comparison in Table 2 is not nearly so favorable. The error bound for 
method A for the example in Table 1 calculated from (5.7) varies from about 10° 
to 10° as x goes from 0 to .5, which is indicative of the fact that error bounds of the 
type we have derived here are generally quite conservative. We note that it is only 
a matter of a little ingenuity to find other examples to make method A appear more 
or less favorable in comparison with other Runge-Kutta methods. 

In conclusion, we emphasize again the main point of this paper which is likely 
to be obscured by giving examples. This is that, if Runge-Kutta methods are to be 
used to start the solution and to change the interval, one is interested only in being 
able to bound the truncation error as well as possible. Thus, that method which 
allows the smallest bound to be put on the error is in this sense best. Therefore, we 
conclude that on a digital computer equations (5.11) and (5.12) should be used 
when a fourth-order Runge-Kutta method is to be used for starting the solution or 
changing the interval. Similarly, equations (4.10) and (4.11) are recommended as 
a third-order system and equation (3.5) as a second-order system. 
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Bounds for Eigenvalues of Tridiagonal 
Symmetric Matrices Computed 
by the LR Method 


By Gene H. Golub 


1. Introduction. In recent years, a number of methods have been proposed for 
finding the eigenvalues of real, symmetric matrices. The methods of Lanczos [8], 
Givens [3], and Householder [5, 14] reduce the original matrix to a tridiagonal 
matrix whose eigenvalues are the same as those of the original matrix. The problem 
reduces then to finding the eigenvalues of a tridiagonal form. Givens has suggested 
the use of Sturm sequences, and others have used Muller’s method [9]. 

In this paper, Rutishauser’s LR algorithm and its variants [10, 11] will be 
considered for tridiagonal symmetric matrices. Henrici [4] has shown that for this 
case the LR algorithm is equivalent to the QD algorithm. It will be shown that 
during the iteration procedure, it is possible to determine bounds on the eigen- 
values. Wilkinson [12] has recently considered the problem of determining rigorous 
error bounds after the eigensystem has been computed. 


2. An Inclusion Theorem. An inclusion theorem is one which exhibits a set 
known to contain at least one eigenvalue. A general discussion of such theorems 
is given by Bauer and Householder [1]. We shall now derive an inclusion theorem 
with the aid of the Lanczos algorithm. 

The Lanczos algorithm for reducing a symmetric matrix A to tridiagonal form 
is as follows: 








Let 
X» = 0, the null vector, 
and 
xX, =X, a non-null arbitrary vector, 
then 
(1) Xeqi = AXe — o4¢Xe — BriXe-1 k=1,2,---,n—1 
with 
, 

@) ~— 
and 

x, Ax 
(3) Ba = —- 


It is not difficult to show [ef. 12] that 


(4) x,’x; =Q for ¢ Aj 
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and hence 
Xn+1 = 0. 


If x lies in a subspace generated by p eigenvectors of A, then 








XpH1 = 0. 
Then since 
x, Axi 1 
Ber oS ee 
X,—1 Xk-1 
(5) as Xe (Xe + ces Xe-1 + Be_oXe-2) 
Xi-1 Xe 
x: x 
k 
- 2 0. 
X,—1 Xk-i 
Now 


(A — ag) Xe = Xess + Brake 
and hence by (4) and (5), 








t_ \2 
= (tt — ov, I)*x, = Xi41 Xi41 + (Xs Xe) . 
Xjy—1 Xe-i 
Thus 
, A ~_ I 2 
(6) x, ( a. )% 6464 = 01 
Xx Xt 
with 


Bo = 0, B, = 0. 


The left hand side of (6) is simply a Rayleigh quotient and as such must lie between 
the smallest and largest eigenvalue of (A — aJ)*. Consequently, 


min (A; — ax)” S ox S max (A; — a)” 
7 i 


where \; , --- , An are the eigenvalues of A. 


THEOREM 1. Let A be a symmetric matrix. Then there is an eigenvalue of A in the 
interval 


(7) u%—-o7SrASatHn. 


If A is a symmetric tridiagonal matrix, then for 


it is not difficult to compute the coefficients a, and 6; . 
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Corouuary 1.1. Let 


aij = a; for t=j 


=bn for |i—j|=1, = min(i, j) 
= 0 otherwise. 
Then for the Lanczos algorithm, 
a= a, 
and 
B= b,” 
and consequently, the interval 
(8) &—-oSASaU+ 


where 
oy = by + bia 
contains at least one eigenvalue. 
Corollary 1.1 is a special case of a known theorem [cf. 2, p. 128]. Note that if the 


intervals are non-overlapping then the bounds given by (8) are smaller than those 
obtained by the Gershgorin theorem. 


If the intervals are non-overlapping, it is possible to obtain improved bounds on 
the eigenvalues by using the bounds of Kohn [7] and Kato [6). 


THEOREM 2. Let 
Ai SASAK OG = 1, 

Then if 

Aj < Aj «. Ajaa 
and if 

Aja < oe < Aju, 
then 

2 2 
(9) mw —-j~*—<rysu4+-—4%-— jf =2,---,n-1. 
Aj+1 — Oe a, — Aja 


Since the coefficients a are Rayleigh quotients, 


\ S min % 
k 


and An 


IV 


Max ay 
A 


In Section 3, the bounds of Theorem 1 will be used. However, the bounds of 
Theorem 2 will be equally applicable. 


3. Application to the LR Algorithm. The LR algorithm proceeds as follows: 
Begin with the given matrix A = Ao, compute a triangular decomposition 
Ao = Loko, 
and then multiply the matrices in reverse order so that 


Ay = Rolo. 








if the 
those 


ds on 
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Now one treats the matrix A; in the same fashion as the matrix A» , and a sequence 
of matrices is obtained by continuing ad infinitum. 
Under certain conditions [10], it is known that A; converges for i+ © to a 














lower triangular matrix with the eigenvalues on the diagonal. 
Let 
(ai? pt O 
1 a 
A; = |. 
(#) 
. al 
O 1 a | 
If the Gauss-Banachiewicz procedure is used to decompose A; , then 
A; = LR; 
with 
| qi >) 
1 
| | 
iF = | ’ 
| 
" 1 P 
1 ei” Dd) 
1 
R; = | 
es 
l ") l j 
where 
qi” es at? 
ab = Bi? /ah? 
gi? = of — of. 
Then since 
Ain = Ry L; 
at = of 449 k=1,---,n—1 
(10) alt) _ 
Bt) = of. gf, 
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If the Lanczos algorithm is performed not with the single vector x; = x but with 
the infinite set of initial vectors xj” = A‘x(i = 0, 1, ---), then it has been shown 
by Henrici [4] that the coefficients a, and 6; may be replaced by a" and @{” and 
that these are the same as those given by (10). Consequently, Theorem 1 may be 
applied to the LR algorithm. 

TueoreM 3. If the LR algorithm is used with the Gauss-Banachiewicz decomposi- 
tion then there is an eigenvalue in the interval 


ay’ — of” SS a4” + of” 
where 
(of”)* = BL” + Bra 
ai? = a? =0. 

Since the coefficients 6;" will become very small as 7 increases, it will also be 
possible to use the bounds of Theorem 2. Theorem 3 enables one to stop the iteration 
procedure when the eigenvalues have attained a predetermined accuracy. 

The LR algorithm converges linearly. Rutishauser [11] has shown it is possible 


to attain cubic convergence if A is positive definite and symmetric. The increase in 
the rate of convergence comes about by working with the matrix 


A; ae yd 
and using the Choleski decomposition. In what follows y; = 0 for all 7 since for 


yi * 0 the eigenvalues are translated, and this in no way changes the argument. 


Let 

















a\® bf? O} 
bo” as” | 
| . . . | 
A; = x 
b\® 
n—l 
lO 1 @.°} 
Then for the Choleski decomposition, 
Aj = R;"R; 
with 
( p\” as? fe } 
p:” 
R; = 
a2, 
e o | 


Dr 





, with 
hown 
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iso be 


ration 


ssible 
ase in 


ce for 
ent. 








TRIDIAGONAL SYMMETRIC MATRICES 443 


where 
(pi?)* = af® 
(dg”)* = (bf°)*/( pt”)? 
(pi? )* = ag? — (df21)”. 
Then since 
Ain = RRZ 
ay” = (pe?)’ + (di?) k= 1,---,n—1 
en” = (ps), 


(O6'*?)* = (di?)*- (PE) 


Note that although a Choleski decomposition of A; was made, it was not necessary 
to compute any square roots. This observation is due to Professor H. Kaiser* of the 
University of Illinois. Since each A ; is a tridiagonal matrix, Corollary 1.1 is applicable. 

THeoreM 4. Jf the LR algorithm is used with the Choleski decomposition then 
there is an eigenvalue in the interval 


a” — oy SS a,” +o 
where 
(on”)* = (be”)? + (04)? 
wi = wi =0. 
4. A Numerical Example. Consider the tridiagonal symmetric matrix A with 
a = 1 k=1,---,5 
b. = $ k=1,---,4. 
When the Lanczos scheme is applied to this matrix, 
a = 1, k=l,---,5 


B. = } k=1,---,4. 


The eigenvalues of A are 


4= 1+ vs = 1.8660254038 
My = 1.5 
As = 1.0 
Ae = 0.5 
M=1l- vA = 0.1339745962. 


* Personal communication. 
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TABLE 1 
F 5 15 
k af of? af” of? 
1 13566009 2.676 X 107 . 13397460 3.255 X 10° 
2 64318606 2.895 X 107 .50011486 7.611.160 
3 1.2211538 4.520 xX 10° 1.0088091 6.813 X 107 
4 1.3750000 4.346 X 107 1.5286215 1.306 X 107 
5 1.6250000 2.602 X 107 1.8284798 Lt xX 
TABLE 2 
k Lower Bounds for , — af” Upper Bounds for \, — a{* 
1 —2.955 X 10° 0 
2 —1.315 XK 10° 1.582 K 10+ 
3 —1.192 k 10° 9.263 X 10° 
4 —9.060 x 10° 3.775 X 10° 
5 0 7.365 X 10° 











The LR algorithm with the Gauss-Banachiewicz decomposition was applied to 
A. The results after five and 15 iterations are given in Table 1. Since the intervals 
still overlapped after five iterations, it was not possible to apply the bounds of 
Theorem 2. However, after fifteen iterations the intervals were non-overlapping, 
and the bounds for \ — a{” are given in Table 2. Note that 


—2.955 X 10° < \, — af” < 0. 


and thus at least seven places of a{” are correct. 
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Polynomial Expansions of Bessel Functions and 
Some Associated Functions 


By Jet Wimp 


1. Introduction. In this paper we first determine representations for the Anger- 
Weber functions J,(az) and E,(az) in series of symmetric Jacobi polynomials. 
(These include Legendre and Chebyshev polynomials as special cases.) If v is an 
integer, these become expansions for the Bessel function of the first kind, since 
J.(az) = Ji (ar). In Section 3, corresponding representations are found for 
(ax)~*J, (ax). Convenient error bounds are obtained for the above expansions. In 
the fourth section we determine the similar type expansions for the Bessel functions 
Y; (ax) and K; (ax). In Section 5, the coefficients of some of our expansions are tabu- 
lated for particularly important values of the various parameters. 


2. Symmetric Jacobi Expansions of Anger-Weber Functions. A function f(z) 
satisfying certain conditions (for these consult [1]) may be expanded in the series 


(2.1) f(z) = >> C0, P(x), -1 S251, a> -—l, 
n=0 


where P,,“~ (x) is called the symmetric Jacobi polynomial of degree n. For our 
present purposes we shall use a definition given in [2]: 





(2.2) 2*nIP,‘* (cz) = (—)"(1 — 2°) *D"[(1 — 2”) **"). 
Also 
(2.3) C.=h, [ @a — 2°)*P,* (x) dz, 
a 2*(n + 1). eves we 
(24) h, = latetDianteat i). ; (v), = “— ’ (v)o = 1. 


Using the representation (2.2) in (2.3) and noticing that. all derivatives of 
(1 — 2”)**” up to and including the (n — 1)st vanish at zx = +1, we integrate 
(2.3) n times by parts to get: 


1 
(2.5) C, = (2"n!h,) / f(z) — 2°)**™ dz. 
1 
Consider the integral definition of the Anger-Weber functions [2, v. 2, p. 35] 


(2.6) J.(ax) + 1E,(axr) = | efe-arsinél dy = f(z). 


0 
When v is an integer, J,(ax) coincides with the Bessel function of the first kind 
J,(ax) [2, v. 2, p. 4]. 
Now differentiate (2.6) n times under the integral sign, substitute the result in 
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(2.5) and interchange the order of integration (which is, of course, permissible) . 
d The inner integral is known [3] and after evaluating it we have 
C, = (—i)*(n + I) a(har”?)* 
(2.7) 
a . —la+(1/2)} 
° [ e“* ees J intetc2(@ sin o) dé. 
ger- 
ials. Use the power series expansion for the Bessel function in (2.7) and integrate term- 
3 an by-term to get 
ince 
for (2.8) C, = (-i)* cos =~ + isin =| A, R,(v, a, a), 
. In 2 2 
ions where 
abu- mae 
a”n! 
A. = 
(2.9) ” ‘ 
r(B+5+1)r(g-s41)n+20+0) 
f(z) 2 2 2 2 
eries and R,, is conveniently described in hypergeometric notation [2, v. 1, p. 182] as 
1 
. Ra(v,a,a) = aFs [+32 41; 
— (2.10) ‘ : 
3 7 v nm ov a 
TS Taeg #5 4+13-$4+5-$] 
Equating real and imaginary parts of (2.6) and (2.1), we get 
(2.11) J.(az) = > AP," (2), —1 S251, 
n=0 
(2.12) E,(az) = > BP, (xz), -1 S21, 
n=0 
where 
Ss of (2.13) A, = A,R, (v, a, a)on (v), 
trate 
(2.14) B, = AR, (v, a, a)¥n(v), 
and 
| (—)*” cos > , n even, 
J} (2.15) on(v) = 
(=O ate > , n_ odd; 
| 2 
, | (—)”” sin =~, nm even, 
k oa , 
ind } (28) ¥.(v) = . 2. 
(n+1) /2 
It in (—) cos —, n_ odd. 


9 ’ 


Equations (2.11) and (2.12) and the expansions in Section 3 may also be de- 
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rived from results in [4]. The present derivation is more satisfactory because it estab- 
lishes a foundation for the work in Section 4. 
When a = —3, 


(2.17) PMO (g) = (4)a(nl)"T.(z), n= 1,2---, 


where 7',(x) is the Chebyshev polynomial of the first kind of degree n. Also for this 
value of a, R, simplifies to the product of two Bessel functions [2, v. 2, p. 11]. With 
a = —}, then (2.11)-—(2.14) become 


(2.18) J.(az) = 2 C,T.(z), -1 S281, 
n=0 

(2.19) E,(ax) = >) D,T,(z), —1 <2 1, 
n=0 

where 

(2.20) C, = & J (ntuy/2 (5) J (n—w/2 (5) onl uv), 

(2.21) Da = nd nt (¢) J (nwa (3) Va(v), 

1l,n=0 
and ¢. “>: 0. 


For integral v we have the expansions 


(2.22) Iulaz) = > 6 Ture ($) | (5) T(z), -l<sr<1 


ll 


, 


n=0 
(2.23) Juri(ax) = 2 » J etn+1 ($) Jaa (3) Tou(zt), -l Sz sl, 
and k = 0,1, 2 --- . Equation (2.22) is known [2, v. 2, p. 100]. 
Since 
(2.24) J.(iz) = e“*? "I, (2), 


where J,(z) is the modified Bessel function of the first kind (2, v. 2, p. 5], we may 
replace a by ia in (2.22) and (2.23) to get expansions for J. (ar) and In4: (ax). 

It is important to note that, although the above expansions are valid only for x 
real and |x| < 1, (2.6) is entire in a and v, and hence a may be chosen arbitrarily 
to yield expansions valid anywhere in the finite complex plane. 

The expansions (2.11), (2.12), (2.18), (2.19), (2.22), and (2.23) are quite 
rapidly convergent, particularly in the Chebyshev cases [5]; consequently the last 
four expansions are eminenily suitable for use on digital computers.* Such series 





* The Bessel functions required to compute the coefficients in our expansions can be 
systematically generated on electronic computers with the aid of techniques discussed in 
(6, 7, 8]. There are numerous tables available for hand calculations. The words “‘accuracy,”’ 
“error,” and “‘convergence”’ in this paper always refer to the properties of the expansion 
when truncated after a finite number of terms. 
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may be truncated and rearranged in powers of z. Clenshaw [9], though, by using the 
recursion formulas satisfied by the Chebyshev polynomials, has formulated a con- 


venient nesting procedure which allows one to utilize such expansions directly. The 
scheme is as follows. Consider 


(2.25) f(z) = >> A,”T,* (2), 0<2<a, 
n=@ a 

(2.26) f(z) = DY An? To. (2) , ~@irsaq, 
n=0 a 

(2.27) f"(z) = D Ae Tro (), —aszsa 
n=0 a 


To evaluate the series (2.25), (2.26), or (2.27), respectively, we construct the follow- 
ing sequences: 





(2.28) b,” = 4 (2) - 2| bots — bs + An”, 
@ 2 
(2.29) b® =| 4 () si 2| b2, — De + An”, 
5 
< 2 
(2.30) bn =| 4 (2) - 2| bats — bate + An”, 
forn = N,N —1,N — 2, --- , 3, 2, 1, 0 with the initial values 
bye: = byae = bya = b¥te = bY}, = OS. = 0. 
f” (x), f° (w), and f® (z) are then given by 
(2.31) f(x) = bo” + d,” [1 —2 (2)], 
a} | 
re 
(2.32) f° (x) = bo” + db,” [1 —2 (2) |. 
(2.33) f(x) ~ [bo _s db,” (2). 


The method is as direct as the ordinary nesting process used to evaluate poly- 
nomials. 
We now derive error estimates for the expansions (2.11) and (2.12) for 
—1<s 2 & 1. Notice that 











(2.34) R,(v,a,a) = 1+0 (2) 
provided all other parameters are fixed, and consequently 
la |" n! 1\| 
(2.35) |An| a+e <a j1+0(2)|, 
r( 5 +1)r( 5) +1) (n+ 20+ 1s 
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and likewise for B, . Also [2, v. 2, p. 206] 

(2.36) max | P,(z) | = (" + *), a> —}. 
—lsz<1 n 


Letey denote the error incurred by taking just N terms of (2.11) or (2.12). Because 
of the rapidity of convergence of the expansions, as shown by (2.35), the (WN + 1)th 
term furnishes us with a convenient error estimate 








a \" Nara 12 1 
| ew | = | | 1 ms 0 = ; 
2.37 a N 
nk gre (Ne 41) (4: *41)ra +1) 
wherea > —3, N>v, -1 S281. 
Among the values of a considered, it follows from (2.37) that the choice a = —3, 


i.e., the Chebyshev case, yields the smallest error term for large N. 


3. Expansions of Bessel Functions of the First Kind of Nonintegral Order. 
Results in the previous section gave symmetric Jacobi polynomial expansions for 
J,(az) and I,(ax) for integral v. When v is nonintegral, these functions are no 
longer entire functions of x, and it is convenient to derive an expansion for the entire 
function 


2.2 
(3.1) T(v + 1)(axv/2) J ,(ax) = oF (o +1;—- =). 


Corresponding expansions for T'(v + 1) (ar/2) “I, (ax) then follow, as before, from 
(2.24). 


Let f(x) in (2.5) be the right-hand side of (3.1). Then we have 


(3.2) Jy(ar) = (ax)” >) APS" (x), —1 Sx <1, 
n=0 
where 
(~)"(ey" 





A, = 
(3,3) QWrl2(Qn + Qa + L)on(m + 4) 0/2) 


‘ 2 
aPa(n +5; u+n-+1, mn +a+3; -¢). 


These equations also follow from a result in [4]. Indeed, using a general expansion 
given there, an alternative formula for (3.3) can be stated. We have 





(34) Fs [-: a; 4 = T(o)(2/2)" (2/ zaSae Pde Ji.-1(2), 


on (—)*29 =P n + 4)(2n +a +4)(Qnt+a+t1). 
e— quite 
(3.5) 





Ea (a/2)*(v + 4), J (a) 
OkiT@ +n+k4+1)° Intk+at(1/2)(@). 
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For the Chebyshev case of (3.2) a = —} and 


(3.6) J,.(ax) = (az)” ¥ C.Pa(2), —-lszzil, 
where 


_— €en( —)"*(a/4)™ 
(37) CG < Qn! T(v +n + 1) 


Notice that when v = —}, (3.3) simplifies. Also, since 





2 
Palm +4; vtn+1,n+1; -¢]. 





(3.8) J-am(ar) = (=) cos (az), 

we infer the expansion 

(3.9) cos (ax) = Eee (x), -l1 S281, 

where 

(3.10) C, = (=) (2n oe pe to 42h J ntarcy2)(a@), 


a formula which can be derived in a number of different ways. 
Using an analysis similar to that of Section 2, we may derive the estimate for the 
error incurred when just N terms of (3.2) are used. 


| | * \a ew | x " g 2@ neta | , + 0 (x) | 
(311) | %! 2er—OONITW +o + 1) (@ + 1) | N/|? 


—-ls2r21, a= -—}, N>v. 





Concerning the optimum choice of a in (3.2), see the discussion surrounding (2.37). 


4. Expansions of Bessel Functions of the Second Kind. The Bessel function and 
modified Bessel function of the second kind are denoted by Y,(z) and K,(z), re- 
spectively, and a treatment of them can be found in [2, v. 2, Ch. VII}. If v is non- 
integral, then 


(4.1) Y,(z) = [sin (vr) {J (z) cos (ur) — J_.(z)}, 
and 
(4.2) K,(z) = 5 [sin (ve) "{I_(2) — 1(2)}, 


so for such values of v expansions for the functions follow directly from the results 
of Section 3. 
If v is an integer, it can be shown that 


43) Wie)@ ae + In (*)| Slee) Modes) < = W,(az), 
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and 


(44) K,(az) = (—)*" E +In (=)| Laz) — = #N,aiaz) +£ WeCiar), 


where 


k-—1 2m—-k 
oF (¢) (k=m—1)! pg 

(4.5) Nin(ar) = Kimmo \ 2 m! 

(0, k == 0, 
and 

- m { @2 we mrt + hm} 

me ~ padies 2 | ) (=) m\(k +m)!" 
In the above y = 0.57721 --- = Euler’s constant and 
(47) he =1$R4+-- +25, ho=d 


We assume the value of log (az/2) is known. Then, since expansions for J; (az) 
and I,(ax) were found in Section 2, and since N;_,(ax) is simply a polynomial in 
1/ (ax), we need expand only the entire part of (4.3), ic., W.(ax), in symmetric 
Jacobi polynomials. 

Using the representation (4.6) as f(z) in formula (2.5), a straight-forward deriva- 
tion gives the series 


(4.8) Wi(az) = > AP," (x), —-1 S251 
n=0 


? 


where 





A, = Ca + ("ln + @ + V(r + @ + 4) 


Qntiatl 


(4.9) > (—)"(—k — 2m). Gy" (hetm + hm) 


& k—n+1 2 m\(k + m)!" 
a+ 2 ntat+l 


2 
We note that the expansion for Y(az) may also be obtained by partially differ- 
entiating (3.2) with respect to v since 


(4.10) Yo(ax) = ayt{ 27efar)} 
du u=0 





A similar procedure yields the expansion for Ky(ar). The Jacobi series for Y; (ar) 
and K,(axr) for k > 0, however, are not so easily obtained in this manner. 
For k = 0 and 1, the Chebyshev cases of (4.3) and (4.4) are 


(4.11) Yo(ar) = = | + In (=)| Jo(ax) + , Ey Ton(x), 0<z 
n=O 


IA 


(412) Y,(ar) = aE + In (*)| Jar) — 22 + D0 Fe T(z), O<2 
1. 2 wax n=O 


IIA 





r) 
ic 


A- 
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(413) Kear) = — E +In (=)| daz) + ¥ Gs Po(2), O<251, 


(4.14) K,(ar) = E + in (=)| I,(az) + = + > A, T(z), e<26 1, 


where 


‘ais _ (3) (I 3 (= :) (n+) hss 


a(n!)? D> (n + 1),(2n + 1), k! 


PP gM 49) et 














(4.16) 
ani(n + 1)! & (n + 2),(2n + 2), k! 
a\” a\* 1 
CO + ie te () = (5) (n + 3), oat 
"(nl & (n+ 1),(2n +1), F! 
a 2n+1 a 2k 3 
(4.18) a Ske (¢) = (3) (n + 3) Tnteti + here) 








ni(n + 1)! i= (n + 2),(2n + 2), k! 


5. Tables. Tables 1 through 3 are based on the Chebyshev polynomial cases of 
the expansions given in the previous sections of this paper. The entries in Tables 1 
and 2 were computed on the UNIVAC 1103-A and those in Table 3 on the IBM 
7090 at ASD. The calculations were designed so that the error incurred in using 
the expansions whose coefficients are tabulated here will not exceed five units in 
the 15th decimal place. Spot checks indicate the error is even less. Because all 
entries are to 16 significant figures, the expansions may be rearranged in powers 
of z with no loss of accuracy. 

The number in parentheses after each entry is the power of ten by which the 
entry is to be multiplied. We have chosen coefficients corresponding to a = 5, but 
the coefficients for other values of a from one through ten are available on request. 

Note that the expansions in this paper are valid not only for —1 S z S 1 but 
for complex z in a region which can be determined by a theorem of Szegé [1, p. 238]. 
More specifically, a Jacobi series representing an entire function converges every- 
where in the finite complex plane. However, the further z lies away from 
—1 <x S 1, the more the accuracy of the expansion deteriorates. This is so because 
P,‘**”(z) for complex x can no longer be bounded by a simple power of n but 
behaves in the following manner [10] 


ping) = MAY) yr (sin 2)" (ond) 


n\_? 2 2 


. cos (NO + 99) {1 +0 (3) 


valid in the z plane cut from —1 to —@ and from 1 to ~. In (5.1), cos @ = z, 


(5.1) 
























































re ee a “a (I—) 92216 LI099 O86T8 "2 sc i eT aR | (gT-) ZIPIP $600 90L48 °Z— rai 
($I—) 10299 989% Zg900°2— (€I—) $8286 S0ZT0 $2802°9 (FI—) 6LZIG Z99FF 8E299°Z (EI—) Z686E 69ZL8 SS86r°z% II 
(ZI—) gge92 16098 ZOEss'*s— (II—) S%Z1¢ #2068 F8I9T'S (ZI—) 89962 6IFIT Sz8El'Z— | (II—) 166ZT S608T T16z8°I— or 
(OI—) P8889 98282 9TOZL'F— (60—) 9LE0F GOLLT S6282'E (OlL—) SOSLb ZSbZE ZO9TF'T | (60—) 16926 09668 62960'T 6 
(80—) 82822 92961 T6E98°Z— (L0—) 8806 FOFEST 80669°T (60—) S9Z8T 68ZIZ TE989°'L— | (80—) SELLL 969LF 8I6hZ°S— 8 
(90—) SZEZT S#986 8Z280°I— (90—) 28098 ZO66I 9ESo6"9 (L0—) ZIIS8 686F% Z0Z0Z'E (90—) OZ68F LIFZL S9TS6'T l 
(SO—) OFZZE SIPEG LIE06°S— (FO—) 980ZL ZSTPO FEBST‘Z (SO—) 680FE LOZSE 6S8Z0°I— (SO—) SL8FF SO6ST 8OFIF'S— 9 
(€0—) 66086 POSOS 60290°T— (S0—) 28829 $1ZZ8 Z6001'S (FO—) Z89Sh LLEIT OL868'Z (S0—) 09266 €F9F9 Z8090°T g 
(ZO—) $2290 S9E8¢ S¢990°Z— (ZO—) 69698 6962Z LZ998'8 (O—) T6868 TOSZO S688z°g— | (ZO—) LE9ZF 9ORET 9EZbE°I— P 
(10—) 69266 S9SSE 862b2°%— (10—) 0F9Z9 SFFL8 SF8T0'6 (ZO—) LL99F SPOTS SFSTO'S (ZO—) ShShL BOOST LL929°6 g 
(00+) S6S1# SZI69 TO88z°2— (00+) 6LFZL I#S6 SEz08"S (10—) S68ZI OLTIO 6I8E9°I— (10—) G9FS9 HOFEZ SzEZz9°Z— z 
(10+) 89689 6900% S6690°'I— | (10+) 600I8 6109 S98h6'T (10—) F2ZE% 18908 SFIIO'T (10—) G8Z8I I#Z89 66999" T— T 
(10+) S@8P9 SZPEL Z2ze88°%— (10+) OFG8E O8S0OL OLSFF'T (10—) 9FFZ9 F90ZE EZL9T'9 (10—) Z988h FEISS 97Z90°2 0 
“H “) “d Sor U 
¢S2z>0 
o=u o=u 
a/t — (g/a) yn I + (ze) {(e/e)u +4} = (zy (g/z)" “MD JI + (2) °T{(8/z) mM + AJ}— = (2)°M 
o o 
o=u o=u 
® w/e — (G/x) pny "I + (2) *¢{(Z/e) up + A}(4/t) = (2)A (g/x)“* "a “ZL + (2)°r{(Z/e)up + 4} (4/%) = (2)°4 
eo eo 
Ss 8914ag ay) Lof szuar2y{o0) 
Z Wavy 
es Pw Z 
(SI—) SZe0¢ SSzzz POSHE'S (SI—) Z609F SSSte Z66hF'T rat 
(¥I—) FP9ZS 96798 BSE2Z°Z (€I—) F619 9BLET SEOs 'Z (FI—) #88E% 100ZE 98ZLe°I— €I—) 006ZZ 22829 88608 °I— Il 
(ZI—) L126 $6692 $9796'T (IT—) SOLS¢ L8S0F SZLPL‘T (ZI—) SI6b6 TI8SZI SP8EI'T ZI—) OL86E $Z90E $8868°6 Or 
(OI—) 20886 62222 SPLIF'T (60—) 9261S $2662 S86FT'T (II—) TOP89 LIZEs8 1z908°2— OI—) 8F890 ZIFSS OFFST'9— 6 
(60—) ITL68 OLI8Z 060ZF 8 (80—) LEZ6T P8068 S899T "9 (60—) LZTS $9086 ITESSe'F 80—) S0FPE OZLE8 806L0°E 8 
(L0—) 20269 96Z9¢ 66080°F (90—) L6S6E 66620 S8FE9°Z (L0—) POSLE O88FE 90226°T— 90—) SE89E T9026 9090Z°I— L 
(SO—) 99LT9 ZE0ZE FHSSTS'T (SO—) $296Z SIFSF 7902'S (90—) $29z9 Sg00S 8zI9¢°9 GO—) 9F69E S9F9E SPEG"S 9 
(FO—) F6EHO FOOTE SPLES'F (€0—) ZI90G L2ZLL SELST'Z ( 60296 10086 $LZ49°T— 99LEE 12992 9%909°L— g 
(80—) L88Zb OT0S6 6ZE90°6 (ZO—) I9T$Z 68026 $¢208°E ( GS6EL OSSIS SLLL8°Z G0I89 SPOTS SZ880'T F 
(1O—) 08628 86699 40608" T (10—) 188 FHE9S FS00S "> ( PESOL ZEI89 SZ96I'°S— E886S 96189 FIESe'6— g 
(00+) S099 99ESS SOTIZ’T (00+) 86999 O8S9T S289z'¢ (T S6II6 96ZET SEZ86°T GSL0Z &2L98 1866'S Zz 
(00+) ISIS¢ SI8I9 O0Szr'9 (10+) 60009 8I€1Z 22992°T (I 800ZZ 09FS9 99F8h* F— 88296 OF860 90Zh6*F— I 
(10+) Zzger 9EzE8 T6S¢9°T (10+) FFFZL SESTH 0&Z80'T (2 ZIZ8E FL86L SZOI8*b— OLSFZ ESZ86 860FE 2 0 
“a “0 “@ “V | u 
2 ¢ 52S g- 
o=u o=u o=u o=u 


(g/2)*"* "qd <4 =(z)'y (¢/2)“*.7%9 s = (ty (G/x) pg 4 =(z)'r = (¢/z)"*"V <4 = (x)/r 


Sarseg 94} 40f 87u9107J20) 
1 @1avy, 





























i Oh 0 SD 


ConmN oO a 














él 
It 
or 


8 (FI—) 8002 Z209F 89F66'T (PI—) 90260 1660 S9T0Z'T 
tet? GITIZ 22986 6LPLL°T (ZI—) LI9ZE STILE 1hez0'I— (PI—) TPOTh OS6IZ OT8E9'°Z (I—) 29616 E9001 29019 I— 
(OI[—) POROF PIOST 9E9ZE"T (II) TO8SE 90FZZ 9O8EZL (ZI—) ZOOSE 9EE99 OFZTZ‘°Z (ZI—) I29%Z% 96ZE6 O8F6Z°T 
(60—) 9Zh6E IPTIOL SGTLI'S (60—) IO8IZ 96188 ZHST'F— (OI—) #1088 19Z98 6ISFE'T (II—) 22686 SF96I S99T9'8— 
(L0—) 98628 SE869 $9820'F (L0-) T&916 02609 819Z6'T (60—) 92699 89928 £8088°8 (60—) 08609 FSFE 68IS9'F 
(GO—) 6P6SE FZIS6 68809'T (90—) 62688 6619% Iz8¢8'9— (L0—) 28880 $0866 ¥2280'°F (L0—) OFIPO 6E862 96986°I— 
a TOLIT 96620 E1898 °F (40—) L6POT 6ESTIO LOSTS'T (SO—) #8h6Z 99TOS F969F'T (90—) #9ZF9 ST899 29z09"9 
ZO—) 96269 29986 99280'T (S0—) #1686 $1906 28z98°e— (FO—) LOO8% O88EE 89T66'S (WO—) 19299 1ZZZzZ 29Z99° I — 
(10—) PPITL OL6PL FRLOL'T (ZO—) ShrOF LHOSD O8PIO'F (S0—) 6662Z $8828 90988" 2 (GO—) TRESS 98z09 e069 'Z 
(00+) 89899 $9662 S8T9L°T (10—) 66020 ZEeez 90z99°Z— (1O—) 69989 L222 669F0'T | (ZO—) IL980 FESS S86EL'Z— 
(10+) SSPZT S2bPE TL060°T (10—) IST9% ZE166 6STZL°9 (10—) $1660 16Z88 P6ZEZ°S | (1O—) SSS00 ZFIOF 6IF8E'T 
(10+) 92906 ZHOTS FEO8E's (10—) 69hIZ Z6SLZ P89ES'Z (00+) 8222 ZESTO ZZS16'S (10—) IIFI6 91Z0E 9Z*Z8"E— 
‘ (10+) IIT LOPPT 4060L°2Z (10—) 908%8 SELSE SEsre*T— (00+) TFPPE 8ZT60 0068'S (10—) SST6S $2061 1Z89Z'T 
2 bs 
= en-)" eie—)"V ee)" (e)"V 
z 
& fat=) 86786 #2900 LP698"°9 (SI—) SOPLL SEEZE S6E6T'F 
a) EE8ZZ PPSLE G9F8Z'9 (€I—) 02820 90299 962L9°E— (FI—) Z80SZ SLZ8Z 62099" 2 (¥I—) SIISt 10008 62F19°b— 
. (IT—) 86068 £8629 622E8 °F (II—) SZO0T 869TF 61269°2 (ZI—) S8ZSIZ 6S9FL 8O8bZ'9 (ZI—) SELL OFSb9 6E669"S 
5 (60—) ZITSZ LOSE9 $8Z80'S (60—) 9IShE 6E9F0 6ZPI9' T— (OI—) #1696 89882 Z61bZ'F (OI—) 06680 S96ET F99TE*Z— 
te) (L0—) #9929 P89ZL L6969"T (80—) Z9886 62928 ST8EL"2 (80—) 90SbZ 9THOZ FORSE’Z (80—) ZhE96 99626 68402" T 
g (90—) 61862 0F999 96999°9 (90—) SO9T8 FI8Zz6 $8006 '°Z— (90—) PIShZ ZISSL SSbHO'T (L0—) 8201¢ Izg08 seee6’b— 
5 t¥0-} 80689 $886Z 28920°2% (SO—) 99020 L696F SE8ZI'8 (GO—) TSI86 86F68 8z909°E (GO—) 8Z0Z9 TPL6h PRLES*T 
- (80—) 0F62% LZ99L SE888"F (G0—) $ZPES FLOPE Sk0Z9'T— (FO—) ZL999 IIPLE Szzge"6 (FO—) 908SF LIZ6L Z098F'S— 
(ZO—) PLES ELISE FHSFI'S (ZO—) SPOTT F2Z69 OF8ZI'Z (ZO—) T2208 26902 988hL°T (€0—) 28862 62460 86968 °S 
an T1ZE9 60862 6ISTO'6 (10—) SIZI8 9ISh6 ZISZ9°'1— (10—) O8@Z1 ShFO0 68P6I'Z (ZO—) 20626 L6FI6 88ELI°S— 
00+) 69ZEL $2268 89Z40'9 (10—) IPLL6 LISTS O&Zbe'S (00+) T8060 TZ60T S8802°T (10—) PI8ES SIFOF $6669 °Z 
(10+) IPShL 22760 ZLZ9T 2% (10—) #8FS0 STSZ9 TZg26°%— (00+) 02629 22908 $000Z°2 (10—) LZ0b F6S8Z IZ618°F— 
(10+) FOSO8 ZEZ6L OFEEL'T (10—) 82928 T98ZE 99TZ0°T— (00+) 986 FITSL L6z8¢"9 (ZO—) SELbT 02629 92066'8 
en)" en—)"V en" en" V 
gcS2z5o- 
o=u o=u 
(9/z)"*Ze"q J = (2) Pe (8/2) Le"V J = (2)*fr-7 
oo oo 
sa1uag ay) 40f 87u9191J00) 
¢ @TavV yy, 
(9I—) 92216 21099 O86TS"Z ary —) ZIT 96600 90z98'%— | Zl 
($I) 10299 9e9¥F zeg00"z— | (SI—) $9286 BOZTO E2802 ] (21) 1¥tr g9v00 soz¥8 ¢— 
a onsen .~ ve £2802 °9 (FI—) 6LZI9 Z99FF 89299'2 (€I—) 26869 69228 S886r°Z 
2 T6098 ZOESS"s (II—) €@ZI¢ F068 FSTOT'S (ZI—) 89962 6IFIT Szser'z— | (II—) I66ZI S608T I16z8°I— 





























(¥I—) 62692 $6906 L1009°Z (FI—) 99086 SZLE ZOT9S*T rat 
(ZI—) LPE8T $8898 £2262 °% (ZI—) SePez 19296 L6Z18°1— (FI—) LEPLG PZ8E6 FISIO'Z (¥I—) S800 S698F O198%° I— II 
(OIT—) 89202 S688Z% FEZOL*T (II—) Lb91Z SFeOb STISZ’6 (ZI—) Z8P8L SLELT STPOL'T (ZI—) 68680 Z68Z% ZI100'T or 
(80—) S8P8E S6ZTS 9L0F0'T (60—) 62Z9TP I8Z86 I8h6z°S— (OI[—) 99828 66922 8ET0Z'T (II—) S0IZ8 FEIS9 889TL'9— 6 
(L0—) 1698% ZS9FS LP9FT'S (L0—) 66169 ZZbPS8 6OZTF'Z (60—) ST&ZE8 LPSLT 092%6'9 (60—) SS6T P6088 ZH6S9'S 8 
(GO—) 8FF9E 216928 60600°2 (90—) #LI6E SEO 9808F°8— (L0—) STO88 8h06F 1Z9zz's (L0—) S8L2Z SETIZ 0Z6L9°T— l 
(PO—) 89829 ShF66 90800°9 (FO—) 92000 80T6I 02802 ‘Zz (GO—) $889% 6ESZT ZOTLTT (90—) L209 99996 99622°¢ 9 
(ZO—) 866b2 98060 SPPZE'T (G0—) 690% 90612 LZ8TO°b— (FO—) Z8E9E 80622 IBLIZ'S (WO—) OTELI ZoSeh ThPLZ°T— ¢ 
(10—) GO8ZE OBL9T 8ZL+0'°Z (ZO—) $8080 86696 FF9L9'F (GO—) SFSS9 966% 0&26E°9 (G0—) ZES06 L9EOF ZZ9ST'Z p 
(00+) 6209% 69998 #8FL0°2 (10—) L8Z%% O9TIT 80296°2— (ZO—) FOLFL IZLE ¥2ZL9'8 (ZO—) 98IZ O1FL9 L6ezE°Z— € 
(10+) 98198 98926 SI8Sz°T (1O—) $26%6 GO00EE FESZ8"9 (10—) 96889 SEST6 sIgse’z (10—) 21662 2699 86988 °T z 
(TO+) PI8S89 06602 FSTHO'F (10—) 22008 80ZSF 28269 'F (00+) $6260 68060 6STIF'S (10—) 26269 FeZIb SSPEs’e— I 
(10+) PZITZ 1602Z LO8TO'S (TO—) #6609 29000 S0TZ0°I— (00+) 16689 Oze0e LO80F's (1O—) FEZ PE9BL 69Z6Z°T 0 
: we)" qwie-)"V qwie)"@ wie)" V u 
& (GI—) P8LZ8 ssoee FISSZ'¢ (QI—) PELES OTZEE ZL8IZ'E (GI—) 99190 L8ZLL OFTHO'T ZI 
y (SI—) 9FL8E 66POS 9L0b8 °F (SI—) TPOTL S996 Ezers'Z— (FI—) S88IZ8 80969 T9986°6 (FI—) 6F68T 009F0 S6r66°S— | II 
(II—) ZOS8Z 28209 OLES2° (II—) TOS9T $8920 +9860°2 (ZI—) LOEST SI6ST ZEZ80°8 (ZI—) FL6ZL 96PO8 $6689" F or 
(60—) ZO068T OT000 LLITH'% (60—) HESS SFSOS 8Z0LZ°I— (OI—) 666% S80Zb 6h8hF'S (OI—) 62682 OF8E8 26096°%— 6 
(L0—) 9818Z 19060 22692°T (80—) S610b S8ZhP ST99T 9 (80—) #286Z ZS0EZ 17866 °Z (80—) TLLZ60 LPL6T 2Esze"t 8 
(90—) GS8ZE S8TFO SELzz'°¢ (90—) ¥998E $98Z0 POEEE'Z— (90—) 62861 ¢89Sz POSTE 'T (L0—) SLZI€2Z 8166Z S6LZT'9— l 
(PO—) $9926 ZE886 8ZFL9'T (GO—) OPOLP 29926 ¥8629°9 (GO—) PL8FI SOShE LEE0S'F (GO—) S9TS¢G 92809 ZIL06'T 9 
(80—) 22689 P6ZPI 66066'¢ (G0—) LhOSS FELPI LEPHE*T— (G0—) TSOZ8 8IZFZ 98PST'T (FO—) 6I08Z LIO9E Z6Sh%'b— $ 
(ZO—) OSL8b SFSSF S00SL°9 (ZO—) 62ZZI 899b0 S9S08*T (ZO—) ZIS9E Se80% LOSZT*Z (g0—) Z908% 69929 PO6rE’9 7 
(10—) 80Z6b 96229 10969" 2 (10—) Z9Z1E L9F6Z 9ZZh° I — (10—) 19808 88%Z6 90ZE9°Z (ZO—) $2229 LEPLb 1E8z0°9— £ 
(00+) 82200 26689 OTS61°S (TO—) 0960¢ $00Z8 8009T “¢ (00+) 1OFZZ YEPET 6OETO'Z (10—) 9FSZF OBEES 88ZI6°Z g 
(00+) SL82Z 9PETL FHLES'T (10—) Leee8 6FSSh 99FbL°S— (OO+) SSS6T OTFFI OfSTE'S (10—) €ZGZ¢ LO8IT SOTL6'b— I 
(10+) LIG8T LPPOb OSFPE'T (ZO—) OTSZ8 86ZTE 660LT'8— (00+) 80890 S080T ZOELF'L (ZO—) 889% 9Z9bZ ZESTZ 2 0 
on-"a on-)"V wo"d wn"V u 














(ponusuoy) ¢ @1av.L, 


& 
+ 








POLYNOMIAL EXPANSIONS 457 


= [n(n + 2a + 1)]"", y = —(1 + 2a)/4. In general, if values of f(x) for com- 
plex z are desired, it is wisest to choose a such that the expansions are interpolatory 
along a suitable ray in the complex z-plane and to stay as close as possible to this 
ray. 

Suppose we have the truncated expansion 


(5.2) f(z) = d, A.Ta() + evir = Ov (Z) teovunr, —LSzrsl, 


and 


(5.3) evi = > A,T, (x) 

n=N+1 
Then ¢y (x) is not generally the Chebyshev approximation of degree N to f(x) in the 
sense of [11], i.e., the polynomial @y (x) of degree N uniquely characterized by the 
fact that in the interval [—1, 1] the number of consecutive points at which the differ- 
ence f(z) — @y(zx) with alternate changes in sign assumes the value 


max | f(z) — y(z) |, 
1sz<1 


is not less than N + 2; but ¢yw(x) may closely approximate #,y(z). How closely, 
of course, depends on the coefficients A, . If A, goes quite rapidly to zeroasn — «, 
then Ay+2 is small compared to Ay; and consequently 


(5.4) ev41 ~ Anus T wii (2) 


and the error curve is practically uniform, i.e., ¢y (x) is nearly y (xz). Such is the case 
in our expansions, and, consequently, we must expect the approximation 4, (zx) 
for moderate values of a to offer a negligible improvement over the Chebyshev 
polynomial expansions derived in this paper and truncated after N + 1 terms. 
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Majority Decision Functions of up to Six 
Variables 


By S. Muroga, I. Toda, and M. Kondo 


1. Introduction. Recently logical elements based essentially on the majority 
decision principle have been widely used in electronic computers. Among these 
elements are parametrons, magnetic cores, transistor-resistor logic, et cetera. 

The logical behavior of such elements can be expressed by a model called a 
“majority decision element” with n Boolean inputs and one Boolean output, whose 
operation can be described in the form of a logical function called a “majority 
decision function’’. 

This paper defines the canonical representative of each equivalence class in the 
classification of the majority decision functions by complementing and permuting 
variables and by complementing the output. Also, a method is proposed to obtain 
all the representatives with their optimum structures, and a table of the repre- 
sentatives of the majority decision functions of up to six variables is provided. 

The reader should be familiar with the content of a previous paper by the 
authors, included as reference [1]. 


2. Majority Decision Functions. A ‘majority decision element” of n variables is 
a logical element with n Booiean inputs, z; , 72, --- , Z, and one Boolean output. 
The output value of the element is 
one for >. wa; = T 
(1)* ry 
zero for do wins i 


IV 


where w; is a prescribed constant real number called a “coupling weight” as- 
sociated with the input z; and T' is also a prescribed constant real number called 
a “threshold.” 

In the case of parametrons or magnetic cores, the coupling weight w; corresponds 
to the number of turns of the winding of the input z; . The threshold T is related 
to the number of turns w, for the constant input by the relation, 


(2) w. = Dwi +1-2T 
where w. > 0 means the constant of one is coupled to the element and w, < 0 
means the constant of zero. 

A set of (n + 1) real numbers (w; , we, --- , Wa ; J’), which specifies the behavior 
of a majority decision element, will be called a “structure” of the element. 

A logical function represented by a single majority decision element will be 
called a “majority decision function.” 


Received September 22, 1961. 
* The term —1 on the right hand side is introduced as a normalizing factor of w,;’s and T’. 
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For example, a majority decision element with the structure (2, 1, 1; 2) repre- 
sents a function zx, + 72273; hence, this function is a majority decision function. In 
contrast, the function x7. + 232% is not a majority decision function since it can 
not be realized by any single majority decision element. 


3. Classification of Majority Decision Functions. Logical functions obtained 
from a given logical function f by the following operations are defined as equivalent 
functions with f: 

(1) Complementation of one or more input variables, 

(2) Permutation among input variables, 

(3) Complementation of f. 

It is a well known fact that the logical functions can be classified into equivalent 
classes by this equivalent relation. Once a structure of a majority decision function 
is given, its equivalent functions can be easily realized in the same element by 
complementing and/or permuting input variables and/or by complementing the 
output. Thus, it is not necessary to determine the whole of the majority decision 
functions; it is sufficient to know the representatives of their equivalence classes. 
It should be noted that this limits the study to a much smailer number of functions. 

In the case of general logical functions, it is difficult to extract systematically 
one representative from each equivalence class, but in the case of majority de- 
cision functions there is a way to define a canonical representative of each equiv- 
alence class from the intrinsic nature of majority decision functions, 

The method of determining the canonical representative is described below. 
Hereafter in this section the majority decision function is assumed to have n non- 
vacuous variables. 

Any majority decision function can be expressed by a polynomial without any 
complemented variable by appropriately complementing one or more variables 
(refer to [1], Section 3). Such a polynomial will be called a “positive polynomial.” 
The way to complement the variables to obtain a positive polynomial from a given 
function is unique if complementing one variable more than once is prohibited. 
Hence we can restrict the possible representatives within positive polynomials. 
This is equivalent to the condition in which the representative should be realized 
by a majority decision element with positive coupling weights. 

All the variables of a majority decision function can be ordered by a relation > 
(refer to [1], Definition 3 and Theorem 1). Therefore, it is always possible for 
variables to be permuted and relabelled so that x: > 12 > --- > 2m holds. This 
permutation can be uniquely determined except in the case of arbitrary permuta- 
tions among some variables such as 2; , 42, -** , Xm for which 7; ~ %2 ~~ +++ ~ Im 
holds. But 2; ~ 22 ~ --- ~ 2m means that the given function is symmetric with 
respect to these variables, and therefore the function is invariant under the per- 
mutations among 2; , 22, --- , m . Thus, the function for which x, > 2 > +--+ > Zn 
holds is unique and can well be adopted as a possible representative. Of course, 
this is equivalent to the condition in which m2 w. 2 --- 2 w, holds for the 
representative majority decision element. Note that as a conclusion from the 
above requirements, we have w, = w, 2 --- = W, > O except w. S 0. 

Only two functions left in each class satisfy both of the conditions just described. 
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If we denote one of them by f, the other is the dual function f* of f. But for 
a majority decision function, either f* > f, or f > f* holds (refer to [1], Corollary 
2). A unique representative of the equivalent class can be determined by requiring 
either of the two inequalities. If we adopt f such that f C f*, this implies w, < 0. 

Thus, it is shown that there is a unique canonical representative in each equiv- 
alent class of majority decision functions which satisfies the following three con- 
ditions: 
Conditions I. 

(1) A positive polynomial, 

(2) m >a °°? Bae, 

(3) f such that f C f*. 

Given a majority decision function, we can now effectively obtain the repre- 
sentative of the equivalent class to which the given function belongs. 


4. A Method to Obtain the Totality of the Representatives of the Majority 
Decision Functions. From Section 5 of [1] it can be determined by linear program- 
ming whether a given function is a majority decision function or not. Therefore, it 
is possible, at least in principle, to obtain the totality of majority decision functions 
by applying the criterion to all of 2” logical functions of n variables. It will, how- 
ever, take an impractically long time to solve 2” linear programming problems 
for large values of n, but the length of time to perform computation will be greatly 
reduced if we can confine the scope of the functions to be tested. 

Accordingly, a method is developed here to obtain a set of logical functions 
which includes all the representatives of majority decision functions and to apply 
the criterion only to those functions in the set. The functions in the set will be 
called “‘candidates” of the representatives. 

Any positive majority decision function can be expressed in the form of Mz, + 
N, where M and N are both positive majority decision functions of (n — 1) var- 
iables, 22, 3, -** , 2n- Therefore, without loss of generality, we can restrict the 
candidates within such functions. This assumes that we have already obtained 
all the majority decision functions of (n — 1) variables; hence the method described 
here is one of the recursive constructions of majority decision functions with re- 
spect to the number of variables. 

Moreover, if we choose as the candidates those functions for which Conditions 
I can be defined, then the set of the candidates will certainly contain the totality 
of the representatives of the majority decision functions of n variables. 

Then the restrictions imposed upon combinations of M and N will be examined. 

Condition (1) will be trivially satisfied, for Mz, + N is positive from its con- 
struction. 

Condition (2) requires that the relation 


(3) 22% 2°: Bt 


must hold for both M and N. Moreover, in order that x; > 22: may hold in Mz, + 
N, it is necessary (Corollary 1 of Reference [1]), that 


(4) mm, 
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where 

m, = M(0, 23, --- , tn) 

m = N(1, 23, --- , tn). 
As the relation > is an ordering relation (Theorem 1 of [1]), the relation 
(5) m1 ZZ *'* ZI 


follows from (3) and (4). 
M and N are majority decision functions satisfying (3), hence the relations 


(6) m 2m, and nm > m 
where 
m, = M(l, 2s, --+ , Zn) 
m2 = N(O, 23, --- , tn) 
hold (Corollary 1 of Reference [1]). From (4) and (6) we have 
(7) MON. 
From (3) in Conditions I, it is necessary that 
(8) f* = N*x, + M*N* Df = Mn +N. 
But as M*N* = M* from (7), (8) reduces to 
(9) M*2>N. 


Thus, we choose as candidates those functions which satisfy the following 
conditions: 


Conditions II 
(1) Both M and N are positive majority decision functions of (n — 1) var- 
iables, t2, %3,--* , Xn. 
(2) For both M and N, x > 43 > -:: > an. 
(3) m 2m. 
(4) M* DN. 


By taking all the combinations of M and N which satisfy Conditions II, we 
can obtain the set of candidates of the representatives of majority decision func- 
tions of n variables. 

M and N must satisfy (1) and (2) of Conditions II. Such functions are either 
canonical representatives of majority decision functions or their dual functions. 
Therefore, once the totality of representatives of majority decision functions of 
(n — 1) variables are obtained, the scope within which functions M and N must 
be taken can be easily determined. In this way we can obtain the totality of the 
representatives of majority decision functions of n variables recursively. 

The next problem is to examine each candidate to determine whether or not 
it is a majority decision function. If so, it is clearly a canonical representative of 
an equivalent class defined in the preceding section. The discrimination of majority 
decision functions from other functions can be accomplished by linear program- 
ming. The details will be found in Section 5 of [1]. 
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5. Majority Decision Functions of up to Six Variables. Following the procedure 
described in Section 4, a program was written for the parametron digital computer 
MUSASINO-I, and all the canonical representatives of the functions of up to six 
variables were obtained. 

The canonical representatives of up to five variables had been obtained by 8. 
Muroga [3] at that time, using a combinatorial method. Both results agreed com- 
pletely. 

The canonical representatives of the functions of up to six variables are shown 
in Table 1. The functions are numbered according to the magnitude of 
V = >°2.1;, which is expected to denote the complexity of functions to some 
extent. Functions are expressed by denoting the variables by means of their sub- 
scripts. For instance, 12 + 13 + 23 stands for the function zz. + 2:22 + Zar. 

In the same entry of the table an optimum structure of the function is shown. 
The optimum structure is one with a minimum number of total turns of windings, 
namely, a structure which minimizes (w; + w. + --- + w, +|w.|) (Section 5 
in [1)). 

To establish the threshold 7, the constant input of zero must be coupled to 
the element with a winding of 27 — V — 1 turns. Dual functions can be realized 
by merely reversing the polarity of the constant input, that is, by coupling the 
constant of one to the same winding. 

The numbers in this table are somewhat different from those shown in [1]. 
This is because f and f* are considered to belong to the same equivalence class in 
this paper and that in Table 1 the numbers of functions of n (nonvacuous) var- 
iables are shown, while the numbers for up to n variables are shown in [1]. 

By computing the number of the members of each equivalent class, the total 
numbers of majority decision functions are obtained and shown in Table 2. 


6. Remarks on the Results. Some remarks are added here concerning the 
representatives of majority decision functions of up to six variables. 

First, it is remarkable that all the candidates proved to be true representatives, 
that is, Conditions II are sufficient for a function of up to six variables to be realized 
by a single majority decision element. 

Second, it is interesting to note that the optimum structures (w;, we, --- , Wa) 
are all integer-valued in spite of the fact that the optimum structure is obtained 
as a solution of a system of inequalities of the form of equation (1). 

A structure of a majority decision function is a solution of a system of 2” linear 
inequalities (Section 5 of Reference [1)}). 


st ie 
Azr=b oof sis pa 


7—-1,2,---,n 
Wi 
(10) Ws 


ors 
Wn 
z 


The third remark concerns the structure of the solution space of these in- 
equalities. It has been noted that for a majority decision function of vp to five 
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TABLE 2 
The Number of Majority Decision Functions 
Number of 
Number of T Number of ae 
Number of Logical Functions a tee a Ma- | Number of Majority | fT" samy 
- of up to n Variables Functions of n .o é of a a Decision 
Variables* a oe « Variables Functions of 
n Variables 
1 4 1 1 2 1 
2 16 2 1 8 0 
3 256 10 3 72 1 
4 65, 536 208 9 1, 536 1 
5 4, 294, 967, 296 | 615, 904 48 86, 080 4 
6 18, 446, 774, ~~ 504 14, 487, 040 | 14 
073, 709, 551, 616 














* These values are obtained from the results in References [4] and [5]. 








TABLE 3 
The Maximum Values of Optimum Parameters of Majority Decision Functions 
n w V = yw T K 
- 
2 1 2 2 | 3 
3 2 4 3 5 
4 3 8 5 9 
5 5 16 9 17 
6 9 33 18 35 














variables the solution space of (10) is a pointed cone. That is, there is a certain 
point x» such that 
(11) Axo = b 


and any solution z of (10) can be written as 


(12) r=t+x Az’ = 0. 


This means the solution space of (10) is a cone with z» as a sole vertex. These 
structures for majority decision functions of six variables were examined and it 
was found that almost all the majority decision functions have solution space of 
a pointed cone but that 15 out of 504 representatives have spaces of non-cone 
structure. These functions are marked with * in Table 1. 

Fourth, some maximum values of the optimum parameters are shown in Table 
3, where V is the sum of coupling weights associated with input variables and K 
is the total number of turns of windings including the constant winding and the 
relation K = 2T — 1 holds. In Table 3, 26 functions have the maximum value 9 
for a weight w and only one function attains the maximum value 33 of V; there 
are 7 functions with maximum K of 35. 
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TECHNICAL NOTES AND SHORT PAPERS 


On the Numerical Solution of a Differential- 
- Difference Equation Arising in 
Analytic Number Theory 


By R. Bellman and B. Kotkin 


Summary. The computational solution of a certain class of differential-difference 
equations requires numerical procedures involving an extremely high degree of 
precision to obtain accurate results over a large range of the independent variable. 
One method of solution uses an iterative procedure which relates the differential- 
difference equation over a large range to a system of ordinary differential equations 
over a limited range. When the characteristic roots of the related system indicate 
borderline stability, it is evident that small perturbations in obtaining successive 
initial values eventually grow out of control as the system increases. 

To investigate this phenomenon, we examine the equation u’(z) = — u(x —1)/z, 
arising in analytic number theory. 


1. Introduction. The function ¥(z, y) equal to the number of integers less than 
or equal to z and free of prime factors greater than y is of obvious interest in number 
theory. It has been investigated by Chowla and Vijayaraghavan, Ramaswami, 
Buchstab, and de Bruijn [1], where references to the other works may be found. It 
has been shown that 


(1.1) lim y“y (y*, y) = y(z) 
yoru 
exists, where y(z) is a function satisfying the interesting functional equation 


(1.2) y (x) = we, a> 6, 
with y(z) = 0, «<0, y(@) =1, OS781. 

The problem of computing the values of y(z) over an initial range, say 
1 S zx S 20, was posed to the authors by M. Hall. Although we possess a method 
described in [2], [3] which reduces this problem to solving successive systems of 
differential equations, the foregoing equation is still of interest, since some useful in- 
formation concerning the accuracy of our technique is obtained. 

Tables of y(x) are given which are of eight significant figure accuracy for 
1 = x S 5 and of two or more significant figure accuracy up to zr = 20. With ad- 
ditional effort, more significant figures could be obtained. 


2. Computational Procedure. As indicated in [2], [3], the single equation of (1.2) 
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TaBLeE 1 
H=2* 
= y(z) y(z + 1) y(z + 2) y(z + 3) y(z + 4) 
1.0000 | 1.00000000} 0.30685282 0.48608400 (10-*) | 0.49109353 (10-*) | 0.354732 (10-*) 
1.0625 | 0.93937538| 0.27701857 0.42592619 (10-*) | 0.42047710 (10-*) | 0.298211 (10-*) 
1.1250 | 0.88221696| 0. 24983249 ©.37274025 (10-') | 0.35958222 (10-*) | 0.250439 (10-*) 
1.1875 | 0.82814974| 0.22504611 0.32575066 (10-1) | 0.30712743 (10°) | 0.210104 (10-*) 
1.2500 | 0.77685644| 020244167 0. 28427227 (10-1) | 0.26199627 (i0-*) | 0.176087 (10-*) 
1.3125 | 0.72806629; 0. 18182739 0.24769805 (10~*) | ©.22521537 (10~*) | 0.147300 (10-*) 
1.3750 | 0.68154627/ 0. 16303349 0.21548879 (16°) | 0.18993579 (10-*) | 0.123315 (10-*) 
1.4375 | 0.63709451) 0. 14590910 0.18716429 (10-') | 0.16141689 (10-*) | 0.103045 (10-*) 
1.5000 | 0.59453490| 0. 13031957 9. 16229603 (10-*) | 0.13701261 (10-*) | 0.860250 (10-*) 
1.5625 | 0.55371290) 0.11614431 0.14050084 (10-*) | 0.11615955 (10-*) | 0.717502 (10-*) 
1.6250 | 0.51449219) 0. 10327493 0. 12143537 (10-') | 0.98366684 (10-*) | 0.597899 (10-*) 
1.6875 | 0.47675186| 0.91613756 (10-) | 0.10479153 (10-') | 0.83206424 (10-*) | 0.497791 (10-*) 
1.7500 | 0.44038422| 0.81072430 (10-') | 0.90292372 (10-*) | 0.70306896 (10-*) | 0.414080 (10-*) 
1.8125 | 0.40529290| 0.71570870 (10-') | 0.77688576 (10-*) | 0.59345175 (10-*) | 0.344148 (10-*) 
1.8750 | 0.37139135| 0.63036275 (10~') | 0.66755423 (10-*) | 0.50041395 (10-*) | 0.285780 (10-*) 
1.9375 | 0.33860152| 0.55402312 (10-') | 0.57290097 (10-*) | 0.42153594 (10-*) | 0.237107 (10-*) 
TABLE 2 
x y (x) 
6 0.196555 (10) 
7 0.87935 (10-*) 
8 0.3640 (10-’) 
9 0.458  (10-*) 
10 0.319 (10-*) 
11 0.284 (10-*) 
12 0.258  (10-*) 
13 0.236 (107%) 
14 0.218 (10-*) 
15 0.202 (10-*) 
16 0.189 (10-5) 
17 0.177 = (10-*) 
18 0.166 (10-*) 
19 0.157 = (10-*) 
20 0.149 (10-8) 





is replaced by the system of ordinary differential equations 


(2.1) 


Equations (2.1) (0) and (2.1) (1) with initial values y(0) = 1, ¥:(0) = 1 are 


(0) 
(1) 


(3) 


(N) 





yo (x) _ 0, 
’ = —yo( x) 
Yi (x) — z+l1 ? 
17) _ Hea (2) 
’ = —Yyn—-i(2) 
yw (x) = z+N 


, 
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integrated numerically to obtain the value y,(1) = y:(0). The three equations 
(2.1)(0), (2.1)(1), (2.1)(2) are integrated numerically to obtain the initial value 
y2(1) = ys(0), etc. Since yi(z) = y(z + DD; j= 0,1,--- ,N, we obtain in 
this way the solution y(r)in OS252N+1. 


3. Numerical Results of Digital Computer Experiments. Using an IBM 7090 
Fortran program with integration subroutines INT and INTM from the IBM 
Share Library D2 RWFINT and a fixed grid size H, our results agreed to eight signifi- 
cant figures up to N = 5 for H = 2” and H = 2~. As N increased, the agreement 
got poorer and at x = 20 there was agreement to only three significant figures with 
the initial value y,(0) = 


4. Stability. Examining the related system of differential equations, we note that 
the characteristic values of the matrix of coefficients are all zero. Consequently, we 
are on the borderline of stability, and progressive loss of accuracy is to be expected 
as N increases. Using finer grids and more precise methods, we could, of course, 
decrease the rate of loss of accuracy. In an earlier computation of the solution of 
differential-difference equations [3] this effect was not present, and more accurate 
results were obtained. 


5. Tables. In this section we present two tables of values. The first presents 
values of y(x) at intervals of 7g accurate to six or more significant figures. The 
second presents subsequent vibiena to the degree of accuracy we possess. Observe 
that as x increases, the number of significant figures decreases. 


The RAND Corporation 
Santa Monica, California 


. N. G. peBruisn, ‘On the number of anaes integers $ : a free of prime factors > y,”’ 


Nedetl. Akad. Wetensch. Proc. Ser. A, v. LIV, n. 1, 1951, p. 1 

2. R. Betiman, “On the computational pot no of differential-difference equations,” 
J. Math. Anal. Appl., v. 2, 1961, p. 108-110. 

3. R. BELLMAN, & B. Kori, On the Computational Solution of a Class of Nonlinear 
Differential-Difference Equations, The RAND Corporation, Paper P-2233, February 1961. 


The Coefficients of the Lemniscate Function 


By L. Carlitz 


Let 9 (u) denote the special Weierstrass @-function that satisfies the differential 
equation 


9p” (u) = 49*(u) — 49 (u). 
Hurwitz [4] put 


P(u) = 





we PE ni 2“ E,, _ 
4 2! 8 6! 4n 4n—2 








foes, 
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He showed that the EZ, satisfy the recurrence 


3 > 4n 
Q) Ee= @—aaeer ay & (Ae - Dn - 4k - 1) (i) Ey En+ 





from which it follows that the Z, are positive rational numbers. 
Hurwitz proved also that 


PY) dain ie? 
u (r+ si)" (4n)! E., 


where the summation is over all complex integers r + si except 0 and 





ont [ = a 
Vi—-z 
The most interesting property proved by Hurwitz for the E, is the following 
analog of the Staudt-Clausen theorem. Let p be a prime of the form 4k + 1 so that 
p = a + 0°; we take a odd and such that 


a=b6b+1 (mod 4). 


Then 4n/(p—1) 
(2) By = G+ 4+ OO, 


where the summation is over all primes p of the form 4k + 1 such that p — 1 | 4n;G, is 
an integer. 

Hurwitz computed the first twelve values of Z, both in the form (2) and also in 
the form 


(3) E, = N,/D, (N,, Dn) = 1, 


where N,, denotes the numerator and D, the denominator of EZ, in reduced form. 

The present writer [2] has discussed arithmetic properties of the coefficients of 
singular elliptic functions. In particular, he showed that if p is a prime of the form 
4k + 3 then 


(4) E,=0 (mod p™) 
when p’ | m, p’ — 1 + 4m and 4m > pr + 1, while 
(5) E, =0 (mod p™™) 


when p’ | m, p’ — 1 | 4m and 4m > pr + 1. In some cases, but not all, (4) and (5) 
predict the correct power of the primes p dividing E,, . However, the formulas give 
no information about the occurrence in the numerator of EZ, about primes of the 
form 4k + 1. An interesting question is whether infinitely many primes 


p=1 (mod 4) 
occur in the numerators of E,, . 


In a later paper [3] the writer obtained congruences (mod 2’) for the E, as well 
as for the coefficients of certain related functions. In particular, he showed that 


(6) % (-" (2) Bete = 0 (mod 2), 


s=0 


provided 2°" |z, r= Oandn > re. 





a 
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In order to get additional information about the numerators N, of Z,., it was 
thought desirable to compute some additional values of N,. The following values 
were computed by R. Carlitz in the Duke University Computing Laboratory, 
making use of the recurrence (1). We remark that the first twelve values are in 


agreement with the results given by Hurwitz. We write N, for the numerator of 
E,, as defined in (3). 


-11°.19°-23*-31°-43?-47°-59-61 -67-71-79-83-859-3137- 
53- 70880471 


In every case the N, has been factored completely. 


It will be noted that N;; = 0 (mod 13) and Ny = 0 (mod 17). This suggests 
the following 


THEOREM. If p is a prime greater than 5, we have 
(7) N,=0 (mod p). 
More generally if 


N, =1 
Ne = 3 
N; = 3'-7 
N, = 3'-7*-11 
N; = 3°-7*-11 
N, = 3'-7*-11"-19 
N, = 3°-7*-117-19-23 
N, = 3”-7*-117-19-23-223 
N, = 3%-7°-11°-19-23-31-61 
Ny = 3"-7°-11°- 197-23 -31 -2381 
Nu = 3°-7°-11'-19*-23-31 
Nw = 3'*-7°-11'-19"-237-31 -43-1162253 
Ny = 3°-7'-11'-19°-23°-31 -43-47 - 13-8887 
Nu = 3-7°-11°-19°-23°-31 -43-47 -G1 -52289 
Ny = 3”-7*-11°-19° *23°-31 -43 -47 -2630966033 
Nig = 37-7°-11°-19*-23-31°-43-47 -59-109-814903 
Nu = 3%-7°-11°-19°-23*-31?-43 -47 -59-17 -80232721 
Ny = 3”-7"-11°-19*-23*-31°-43 -47 -59-67 -48316510111193 
Ny = 37-7"-11°-19*-23°-31°-43 -47 -59-67-71 -3469- 1330177 
No = 3°-7"-11'-19*-23*-31°-43 -47 -59-67 -71 -503 - 1389248989981 
Na = 3"-7"-117-19*-23°-31°?-43-47-59-67-71 -79-221430324996967 
Ne = 3"-7"-11°-19*-23°-317-43-47 -59-67 -71-79-83 -6195097 - 111333763 
Nx = 3°-7"%-11°-19*-23*-317-43?-47 -59-67 -71-79-83 -474230496549203 
3".7" 
77 


(8) p\|m, p— 144m, (—1)"2""" #1 (mod p) 
then 
(9) N,=0 (mod p’). 


Proof. lf we put 


4m 
U zx ‘ \—1/2 
Eg EE = (@(u))-”) 
o(u) = Be (4m)! o(u) (Pu) 
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then, as Hurwitz showed, 
(10) Bim = (1+ 2)*{ (1 + 1) — QE... 
The writer [1, Theorem 9] has proved that if p’|m, p — 1 + 4m, then 


Bim =O (mod 7’). 
Thus (10) yields 


(11) {1 +14)" —2}E,=0 (mod p’). 
Since p — 1 + 4m, the denominator of £,, is not divisible by p, so that (i1) implies 
(12) {ai +i)" —2}N,.=0 (mod 7’). 


In the next place, since 
(1 + 4) = (—4)", 
it is evident that 
(1+ 4)°"—2=0 (mod p) 
if and only if 


(13) (—1)"2"""" =1 (mod p). 
If, therefore, (13) is not satisfied, it is clear from (12) that 
Nn =O (mod p’). 
Finally we note that if p = m, then (13) is not satisfied and (7) follows at once. 


Duke University 
Durham, North Carolina 
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2. L. Carurtz, ‘“The coefficients of singular elliptic functions,’’ Math. Ann., v. 127, 1954, 
p. 162-169. 


3. L. Caruitz, ‘Some arithmetic properties of the lemniscate coefficients,’’ Math. Nachr., 
v. 22, 1960, p. 237-249. 

4. A. Hurwitz, “Entwickelungscoeffizienten der lemniscatischen Funktionen,’’ Math. 
Ann., v. 51, 1899, p. 196-226 (Mathematische Werke, Basel, 1953, v. 2, p. 342-373). 


All Factors g < 10° in All Mersenne Numbers 
2° —- 1, p Prime < 10* 


By H. Riesel 


During the year 1960 the author made additional investigations respecting 
factors of Mersenne numbers M, = 2” — 1, where = is a prime. The author has 
earlier examined the least factor g of M,, for p < 10°, if gq < 10-2” ~ 10’. This 
inquiry was secondary to an effort to ascertain Mersenne primes [1]. 

The present examination, which resulted in 255 new factors of the numbers M, , 
was done in the following manner: All primes s < 10°, for which s = 1 (mod 2p) 
and s = 1 + (mod 8) were tested as factors of M, ter 1200 < p < 10000. The 
range p < 1200 has been examined by Brillhart and Johnson [2]. 
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To save computing time the following scheme was used. If it takes long to decide 
whether a number s of the kind in question is a prime, but an examination of 
whether 2? = 1 (mod s) can be accomplished in a short time, it may be convenient 
to test unnecessarily some s-values which are not primes. Such values can be 
obtained by using such s-values as do not contain any factor St, where ¢ is given. 
If +/s < t < s only primes remain. Generally, if t < +/s the fraction 


: 1 - 

y= v(t) -11 (1-4), p prime 
p=3 Pp 

of all numbers remaining after this test. Preparing the program for use on the 

BESK, it was found by test runs that the computing time was minimal with the 

following t-values: 


= 40 for 1200 < p < 2000 
70 for 2000 < p < 3000 
90 for 3000 < p < 4000 
120 for 4000 < p < 5000 
130 for 5000 < p < 6000 
170 for 6000 < p < 7000 
= 200 for 7000 < p < 10000 


The running time on the BESK was about 20 minutes for each p. 
All Mersenne numbers M, with p < 10000 have now been examined with respect 
to all factors g < G = 10°; for some M, this limit G has been considerably raised. 


A summary is given below of what in the author’s opinion has so far been accom- 
lished in this field: 


Ce ee 
ll 


ps 97 completely factored 

p = 101 G=2* 

103 < p < 149 G = 2" 

p = 151 completely factored 

157 S p S 257 G = 2" 

263 < p < 389 G = 2” 

p = 397 G = 2 

401 < p < 1021 G=2 

1031 <ps51193 G=2" 

12001 < p< 2999 G= 10° 

3001 < p < 3593 G = 10°, if any factor is known 
3001 < p S 3593 G = 38400p + 1, if no factor is known 
3607 < p< 9973 G= 10° 


All hitherto known factors are to be found in [1], [2] and this paper. Furthermore, 
the author has calculated the two largest factors of Ms; and Ms; ; these are 


13842607235828485645766393 and 7289088383388253664437433 


respectively. 
As has earlier been pointed out [3], two typographical errors occurred in [1]; the 
factors belonging to Mogg and Msus should be 7158119 and 643217, respectively. 
Some of the factors referred to below have been calculated by Edgar Karst [4], 








480 H. RIESEL 


[5], namely, the 24 factors designated by an asterisk(*). Furthermore, all factors 
stated in [1] and [2] as well as in this paper have been checked with a special pro- 
gram and have been found correct, except for the typographical errors mentioned 
above. 

Lastly, the author desires to express his gratitude to the Swedish Board for 
Computing Machinery, which has graciously provided machine time on the BESK 
for this study. 


Stockholm-Vallingby 
Sweden 


EDITOR’S NOTE: During the period in which this paper was in process of 
publication A. Hurwitz [6] has listed 92 Mersenne composites, corresponding to 
primes p in the interval 330 < p < 5000, for which no factor is known at present. 
He also lists Mas; as such a composite number, and he announces the primality of 
Mas; and My;. 8S. Kravitz [7] has extended the writer’s earlier table [1] to the 
range 10,000 < p < 15,000. 
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2. J. Bruutwart & G. . JOHNSON, “On the Factors of Certain Mersenne Numbers,”’ 
Math. Comp. v. 14, 1960, 365-360 
. JL. SELFRIDGE, TE 271, MTAC, v. 13, 1959, p 
t E. Karst, “New factors of Mersenne numbers, * Math. Comp., v. 15, 1961, p. 51. 
5. E. Karst “Faktorenzerlegung Mersennescher Zahlen Mittels Programmgesteuerter 
nee Aum ” Numer. Math., v. 3, 1961, p. 79-86. 
. ALEXANDER Hurwitz, “New Mersenne primes,”’ Math. Comp., v. 16, 1962, p. 249-251. 
. Smpney Kravitz, “Divisors of Mersenne numbers 10,000 < p < 15, 060, om” Math. Comp., 
v. 1K 1961, p. 292-293. 





P Factors of M, p Factors of M, 
1201 1967239 -8510287 2069 1816583 

1223 31799 2081 3308791 -3874823 -7920287 
1229 46703 2099 22858111 

1231 5684759 2113 52030513 

1249 52358081 2131 17359127 

1321 1394977 -4848071 2251 778847* 

1361 3397057 2309 61636447 

1367 65617 -232391 -2561759 2339 299393 -1445503 
1399 28875361 2389 172009 -267569 
1433 53101249 -86576129 2393 10758929 

1439 46049 -172681 2399 5201033 

1543 472159 2411 1393559 -4185497 
1583 189961 -8589359 -43817441 2579 19868617 

1627 25722871 2593 62233 

1637 81679753 2609 7404343 

1663 6013409 2617 2213983 

1667 2493833 - 10975529 2621 12895321 

1723 56421359 2657 148793 -318841 
1741 4230631 2663 15743657 

1789 254039 2677 465799 

1861 2735671 2687 1918519 -2434423 
1879 13697911 2689 90092257 

1901 12147391 2699 307687 -1187561 
1931 50207 2711 76954447 

1979 32392273 2731 93968249 

2063 53639 -165041 2749 45737863 
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P 
2843 


2917 


3037 
3067 


3119 
3121 
3167 
3181 
3191 
3253 
3257 
3299 
3329 
3391 
3433 
3529 
3557 
3593 
3701 
3719 
3761 
3767 
3779 
3793 
3863 
3923 
4003 
4021 


4127 
4139 
4211 
4373 
4507 
4597 
4759 
4903 
4919 
4931 
4943 
4957 
5003 


5087 
5119 
5189 
5197 
5227 
5231 
5273 
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Factors of M, 


187639 

60457879 

22869281 

1374833 -68747527 
95278231 

145777* 

5565031 * 
22063999* 
15914447 * .68886553 
230807 * - 14222641 * 
31509617* 
12237289* 

127241* 
40895857 * 
46452841 * 
4032167 * 
19873177* 

665801 * -1005359* -26225863* 
1519169* 

5952823* -12688369* 
5335849 * -13523129* 
2952311* 

6086543 * 
318287 - 1561823 
9490889 

55655279 
50455199 

2086009 -6084191 
91033 

1120271 

2236111 

16756559 

10800407 

262337 

2080009 

16820897 

92271433 -94528529 
104953 

21381209 

73553 -93769607 
161807 

13973551 

2872697 

33264527 

128519 

3003943 

1050631 

36425449 

1678711 

82978991 

3362473 

20195543 

47126633 

41849 -470791 
38134337 


P 
5281 
5303 
5381 
5417 
5431 
5507 
5521 
5563 
5569 
5639 
5717 
5741 
5821 
5827 
5867 
5903 
6067 
6121 
6199 
6271 
6287 
6311 
6317 


6373 
6427 
6481 
6491 
6529 


6761 
6841 
6871 
6947 
6977 
6983 
7079 
7159 
7187 
7211 
7253 
7349 
7451 
7457 
7499 
7517 
7547 
7673 
7717 
7841 
7883 
7901 
8087 
8111 
8123 


Factors of M, 


37813511 
2447327 

45737 
48224401 
88805177 
91658711 
16239857 
15489473 -42182839 
3118439 
81237913 
743881 
13959247 
5180489 -57186553 
68423863 
4548241 
2172727 
140207 
10951609 
2799793 -38678609 
69440719 
1201337 
26777041 
26340857 
64962137 -77973719 
39027281 
333457 

209311 
13630817 
56633 - 1486591 
26058761 
48871601 
57689 -475927 
52613263 
42197959 
95506919 
19969847 
99121783 
12012167 
16980751 
13918823 
19446841 
6837353 
441449 

442457 -679487 
96671999 
41917649 
25798649 
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Pp Factors of M, p Factors of M, 
8167 76835137 9049 28721527 -28938703 
8171 9412993 9059 30293297 
8209 14759783 9109 49625833 
8221 9667897 - 18480809 9127 8707159 
8269 19630607 9137 2704553 
8273 28062017 -62014409 9161 86901247 
8287 36877151 9199 53354201 
8377 134033 -787439 -2596871 9221 91841161 
8429 455167 -927191 9283 29352847 -34031479 -41532143 
8467 6655063 9337 2838449 -2405633 
8539 13662401 9403 5735831 
8563 32402393 9479 48532481 
8573 12345121 9601 3513967 -16974569 - 17256487 
8623 80504329 9643 12362327 
8699 43790767 9743 34626623 
8737 6640121 F 9817 20556799 
8741 5926399 ; 9829 14075129 
8849 52368383 — 9851 3723679 
8933 36232249 9859 1656313 
8969 13345873 9883 10436449 
9029 25913231 9973 7419913 -10591327 - 19367567 


Certain Properties of Pyramidal and 
Figurate Numbers 


By M. Wunderlich 


It is well known that despite some extensive computation [1], the only two 
known solutions to the Diophantine equation 


(1) a+b+c=3 
area = b= c = 1, anda = b = 4,c = —5. Professor Aubrey Kempner noted at a 


number theory seminar at the University of Colorado that these solutions also 
satisfy the equation 


(2) a+ +cec=a+b+e. 


Therefore, it is of interest whether or not (2) has solutions other than these two 
and if so, how many. Since there are so few solutions known to (1), it seemed 
reasonable to conjecture that there would be only finitely many solutions to (2). 
If we change the sign of the third variable and divide through by six, we see 
that (2) is equivalent to 
a@—-a -b c-c 


ey Wg 





or 
(a — 1)(a)(a + 1) 4 (b — 1)(b)(b + 1) _ (ce — 1)(e)(e + 1) 
6 6 6 ; 


Received March 13, 1962. The research for this paper was supported in part by the National 
Science Foundation. 
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Numbers of the form 


had been studied by the ancients and had been given the name “Pyramidal” because 


of a curious geometric property which they possess [2]. It is also clear from (3) that 
we can define the pyramidal numbers, P, , as follows: 


AES ictteg | n=l. 


In view of the binomial theorem, they are the numbers in the fourth diagonal row 
of Pascal’s triangle, a fact which makes the numerical computation of these numbers 
very easy. 

Therefore, the conjecture mentioned concerning (2) can be restated as follows: 
There is only a finite number of solutions to 


(4) P,+ P, = P.. 


The purpose of this paper is to indicate how the reasonableness of this conjecture 
was tested and how an examination of a table of solutions of (4) led to some inter- 
esting theoretical results. In particular, S. Chowla was able to prove that there are 
infinitely many solutions to (4) and hence to (2) [3]. Also, S. Segal could prove 
that the only solution to 





a=2 


2P, = P, 
is where x = 3 and y = 4 [4}. 

The computation of the first 88 solutions of (4) was done with the help of V. 
Keiser on the C.D.C. 1604 digital computer at the National Bureau of Standards 
laboratories in Boulder, Colorado. The method used was systematically to deter- 
mine for every pair of pyramidal numbers P, , P,,, m < 13,000, n < m, whether 
or not P,, — P, was again pyramidal. If we let M = P,, — P, it is necessary to 
determine whether 


M=P,, v= 1,2,---n—1. 


To do this, the following “hunting” procedure was employed: All the numbers 
P, , P:, --- P, were stored in the machine in ascending order of magnitude. M was 
first compared with Pj,,;2;. ({m] as usual indicates the integral part of n.) If M = 
Pinj , @ solution to (4) was found. If M < Piq,; , a solution can only exist for 
v = 1,2,--- , [n/2] — 1. If M > Pia , a solution can only exist for v = [n/2] + 
1, [n/2] + 2, ---., n. In either case the number of values which v can assume is 
roughly halved. If this procedure is repeated [log m/log 2] + 1 times, any solution 
if it exists will be found. For example, if m = 15,000, the process need only be 
repeated 14 times. 

Table 1 lists the first 88 solutions found using this program. It required approxi- 
mately 6 hours of machine time. 

It is interesting to note that although S. Chowla proved that there exist infinitely 
many solutions to (4), he by no means justified the great number of solutions that 
were found. By imposing extra relations on (2) he was able to reduce it to a Pellian 
equation which has infinitely many solutions that also satisfy (2). Therefore, he 
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TABLE 1 

x y Zz = y z 

3 3 4 1351 1478 1786 

8 14 15 798 1818 1868 
20 54 55 438 2164 2170 
30 55 58 1146 2072 2183 
39 70 74 1139 2115 2220 
61 102 199 1609 1941 2256 
84 90 110 1105 2303 2385 
34 118 119 853 2417 2452 
48 138 140 1103 2514 2583 
119 154 175 1484 2584 2738 
187 201 245 1089 2773 2828 
100 290 294 834 2958 2980 
327 336 418 528 3138 3143 
149 429 435 1775 2954 3154 
252 424 452 1484 3094 3204 
248 450 474 2478 2726 3286 
362 415 492 2099 3211 3486 
219 515 528 729 3595 3605 
136 532 535 2200 3660 3908 
424 448 550 742 4415 4422 
314 527 562 2116 4580 4726 
434 495 588 2948 4408 4810 
399 588 644 3138 4630 5068 
324 663 688 2912 4838 5167 
272 688 702 868 6034 6040 
304 695 714 2252 6390 6482 
349 713 740 5338 5608 6900 
532 643 747 3570 7154 7439 
424 705 753 1271 7554 7566 
378 790 818 6152 6586 8034 
608 754 868 1160 8070 8078 
230 903 908 5300 7284 8120 
489 869 918 5630 7105 8129 
775 950 1098 6340 6788 8280 
703 1064 1158 4115 8034 8379 
878 1044 1220 4015 8910 9174 
968 1001 1241 7104 7847 9442 
922 1286 1428 7062 8094 9592 
290 1430 1434 2951 10184 10266 
367 1436 1444 1328 10568 10575 
855 1343 1450 7842 10168 11532 
504 1629 1645 7294 10618 11660 
897 1621 1708 8274 10149 11725 
750 1690 1738 9050 11100 12824 








found infinitely many solutions of a very special type, none of which, incidentally, 
appear in the table.* Two unsolved problems are to find a parametric representation 





* (Added in publication) In the March 1961 issue of Elemente Der Mathematik, W. 
Sierpifiski has also shown that there are infinitely many such solutions. His proof, however, 
does yield two of the solutions in the table, namely z = 8, y = 14,2 = 15 and z = 2912, 


y = 4838, 2 = 5167. 
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which will give all the solutions to (2) and to find an asymptotic density function 
analogous to the prime number theorem. This latter problem can be more explicitly 
stated as follows: Let ¢(z) represent the number of integers n < z such that there 
exist positive integers a, b such that 


P.+ P= P,. 
Does there exist a continuous function g(z) such that ¢(z) ~ g(x)? That is, 
_ o(2z) 
lm —— = 1 
zo g(x) 


The concept of pyramidal number can be generalized by defining the rth figurate 
number of order n to be the binomial coefficient 


(5) het tS. 


n 


In this notation the pyramidal number P, is f;,. . Work is now in progress to com- 
pute possible solutions to 


Sas + Jas ” Jus 


where n = 4, 5, and 6. One might be led to believe that there are only finitely many 
solutions to (5) for n = 4 for the following reason: Whereas for n = 3 the equation 
was reduced to a Pellian equation, if nm = 4 & similar reduction may result in a set 
of cubic Diophantine equations to which the Roth theorem may apply. If the reduc- 
tion could be effected without imposing any further restrictive relations the con- 
jecture would be proved. Preliminary results show that the only two solutions to 
(5) for n = 4 and z < 5264 are z = 4,y = 4,2 = 5; andz = 129, y = 187, 

= 197. It is interesting to note that Paul Erdos had once conjectured in a letter 
to Mr. Chowla that the only solution of 


2fn.z ss Saw 


isz = nand y = n + 1. As we have seen, S. Segal has affirmatively confirmed this 
conjecture forn = 3. 


It was further noted upon examining a decimal print-out of the first 25,000 
pyramidal numbers that the last digit repeated itself in a cycle of 20, i.e., 


(6) If z=y (mod20), then f3,. = fs, (mod 10). 


This observation led to the following generalized result: 
Turorem. If k = p,"'p.™ --- p,"*, and n is a positive integer then for j = 1, 


2,--- qlet 
Ell +E 
Pi Pi Pi’ 


where p;*' is the largest power of p; S n. i.e. a; = [log n/log p,). Finally let t = pi” 
pi? --+ po. If x and y are positive integers such that x = y (mod tk), then 


Jax = | (mod k) 
(Note that (6) is a special case of the theorem where n = 3 and k = 10.) 
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Proof. For i = 1, 2, ---n, (© +n —i) = (y+n-— i) (modtk) so that 
(c-+n—1)(a4+n—2) --- (x) =(y+tn—1)\(y+n—2)---(y) (modék). 
Also from the definition of ¢, (n!, tk) = t since ¢ is the product of highest powers 
of 71 , P2, *** Pq Which are contained in n!. So, 


- oes) -:: (z) _ (y+n—1) --- (y) 


n! n! 





” fny(mod k) 


which proves the theorem. 

It is not asserted in the theorem that tk is the smallest period. In fact, easy 
examples show that in many cases a value for ¢ can be found which is strictly smaller 
than the one specified in the theorem. According to Mathematical Reviews (v. 20, 
1959, Review no. 1653), the smallest period has evidently been found by 8. Zabek 
[5] to be tk where 

t= pr” po --+ pa’®. 
This information may be quite useful in numerically searching for solutions to (5) 
since these congruences limit the number of solutions that could possibly exist, 
thereby reducing the amount of machine time needed for the search. 


The University of Colorado 
Boulder, Colorado 
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Note on Osculatory Rational Interpolationt 
By Herbert E. Salzer 


Abstract. In n-point osculatory interpolation of order r; — 1 at points z;, 
t= 1, 2, --- , m, by a rational expression N(z)/D(x), where N(x) and D(z) are 
polynomials >> a,’ and >. b,’, we use the lemma that the system (1) 
{N(2:)/D(a:)}™ = f (2:),m = 0,1, --- ,r; — 1, is equivalent to (2) N“”(x;) = 
{f(2i)D(x;)}, m = 0,1, --- , rs — 1, D(a) ¥ 0. This equivalence does not re- 
quire N(x) or D(x) to be a polynomial or even a linear combination of given func- 
tions. The lemma implies that (1), superficially non-linear in a; and b; , being the 
same as (2), is actually linear. For the n-point interpolation problem, the linear 
system, of order >+7-: r;, which might be large, is replaceable by separate linear 

Received April 20, 1961. 
t Much of the material in the present note is contained in an entirely independent (still 


unpublished) study by Henry C. Thacher, Jr., who was kind enough to send the writer a copy 
of his preliminary draft. 
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systems of orders r; (or even 7; + ris. + --- + ri; When conveniently ernall) by 
applying the lemma to the continued fraction (3) N(z)/D(z) = a» + o—*| 




















| daa 
_ =— I—2 _ ~ _ 
ad Mla 4 2 | il,2 tly... 4% w|, 2 tly. 
| 4,2 | Q,r,—1 | a2,0 | 2,1 | da,r2—1 | 43,0 
wel l = me 
ye eee) FS a | + -- $e | . In (3), which has the property (proven 
| an,o | Qn, | Gn,r,—1 


in two ways) that the determination of a;.,, is independent of all a’s that follow, we 
find a;,, stepwise, but several at a time (instead of singly which is more tedious), 
retrieving them readily from the solutions of those lower-order linear systems. 


1. Introduction and General Lemma. In n-point osculatory rational interpola- 
tion for a given function f(x) by N(z)/D(x), where N(x) = >> aw’ and D(z) = 
> bs’, we have to find the coefficients a; and b; from the following d = 5°?) r; 
conditions: 


; = ) = 0,1,---,r—1, 
(1) d™{N(x)/D(x)}/dx™ | a2, = f(z), 


4#=1,2,---,n. 
The N(x) and D(z) are usually taken to be of nearly equal degree, i.e., for odd d, 
both of degree [d/2], and for even d, of degree [d/2] and [d/2] — 1 respectively. 

One answer to (1) is provided by Thiele’s reciprocal difference formula for 
confluent arguments [1]. But that requires the build-up and tabulation of a recipro- 
cal difference scheme involving confluent forms, which might be too cumbersome to 
handle for a large total number of conditions. The present approach considers a more 
direct solution of (1), and avoids confluent reciprocal differences. 

Whenever in (1) some r; > 1, the determination of a; and b; after differentiation 
of the left member, appears offhand to involve the solution of a non-linear system of 
equations. In reality, (1) may always be solved by an equivalent linear system. 
Before citing the general lemma which establishes this equivalence, it is instructive 
to verify the first few cases. Let N and D denote any functions of x, not necessarily 
linear combinations, the subscript 7 denote the argument x; , and D; ~ 0. 

Ordinary rational interpolation, (N/D); = fi , is, of course, equivalent to N; = 
(Df); - 


First-order osculatory interpolation, 
(2) (N/D)i=fi,  (N/D)i = fi, 


may be expressed as N; = (Df); and N//D; — N.D//D? = f/. In this last equa- 
tion, replace N;/D; by f; and multiply by D; , so that (2) is replaced by the equiva- 
lent 
(2’) Ni=(Df)i, Ni = (Df). 
Whenever WN and D are linear combinations of specified functions, the system (2’) 
is linear in the coefficients. 

Similarly, for second-order osculatory interpolation, namely, 


(3) (N/D);: = fi, (N/D)i = fi; (N/D).” = fi’ 


= <> 


carrying out the differentiation in the last equation of (3), in view of the equivalence 
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of the first two equations in (3) to (2’), we find (3) to be equivalent to the system 
(3’) Ni = (Df)i, > Néi= (Df)i,  Ni& = (Df). 


The foregoing equivalence relations suggest this general lemma for any N and D, 
as long as D; ~ 0: The system 


(4) (N/D)$” = f$”, m = 0, 1,---,7, 
j8 equivalent to 
(4’) NS” = (Df), m=0,1,-+-y7. 


The proof is immediate by induction. Assume the equivalence of (4) and (4’) up 
to order r — 1. Then first assume (4’), and apply Leibnitz’s theorem to both (4’) 
just for m = rand to N$” = {(N/D)D}§”. In the former, replace f{” by (N/D)$”, 
m = 0,1, --- ,r — 1. Comparison of terms shows that (N/D){” = f$”, so that (4’) 
implies all of (4). Conversely, starting from (4) and applying it to the Leibnitz 
formula for N$” = {(N/D)D}§$°", we obtain the equation in (4’) for m = r, which 
establishes all of (4’). 

A discrete case analogue of this lemma (i.e., for finite differences or divided 
differences, instead of derivatives) is established even more readily. Thus the equal- 
ity of the differences of N/D and f at xz; implies that (N/D),; = f; at every 2; of the 
total > r; points (assume now that every D; ~ 0), from which follows the 
equality of N; and (Df); which, in turn, implies the equality of the differences of 
N and Df at every zx; . The lemma above is seen to be a limiting confluent case of 
the finite difference analogue. 

Application of the lemma to (1) reduces the general problem of osculatory 
rational interpolation to the solution of just a linear system for the a; and b; , which 
is thus not harder than the problem of osculatory polynomial interpolation. It is 
apparent from the inductive process in the proof of the lemma that (unlike the 
polynomial case) we could not expect to obtain a linear system from any modifica- 
tion of (1) where m either fails to start at 0 or to run consecutively, even at just a 
single point 2; . 

By this lemma, even when N(x) and D(z) are not polynomials, but linear com- 
binations of preassigned known functions ¥,(x) in place of x’, we may replace the 
system (1) by the equivalent linear system 


, (m) = 0,1,---,n-1, 
(1’) N™ (zi) = d™{D(x)f(x)} /dx™ | ome; 5 
ye ee 


where every D(z;) ¥ 0. 


2. Application of Lemma to Continued Fraction Interpolation. From now on 
we shall consider just the most important case when N(x) and D(z) are poly- 
nomials. The lemma may be applied to solve (1) or (1’) by a number of linear 
systems, each of much lower order than that of (1’), namely, >.7.1 r;, which might 
be inconveniently large, by the method described below. 

We may express N(x)/D(z) in (1) as the continued fraction 





ieee 











(5) 
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N(x)/D(z) » ais ¢ Sel peme) vu 422" 4 eel 


| 4,1 | 4,2 

















| Q1,r;-1 | a2.0 
(5) Odie --> + tect dee yas 
| dea | @e,rs—a | ds,0 
4 2a teal T= tel ore 47 — al 
| @n,o | Qn. | Qnrn—1 


Now (5) has the useful property that each coefficient a;,,, which may be found from 
(1), involves just the preceding coefficients a;, and is completely independent of 
every following a;, . This property may be seen from the limiting confluent form of 
Thiele’s continued fraction formula in terms of reciprocal differences [1]. Another 
direct way to establish this property without employing the Thiele reciprocal dif- 
ference scheme, is by induction upon (5). Thus, at any z; in (1) and (5), the in- 
terpolating equation for the next higher order derivative, say the mth, may be shown 
to involve a;,, a8 the only new quantity, provided that the same is true for all 
coefficients preceding a;,m , likewise down to a;,-,1, then a;4:.9 is found from the 
value of f(zi4:), and so on down to a,,,,-, . By inspection, this property holds for 
the first few coefficients a;,o , a1,1 , which suffices to complete the induction. 

To apply this property, we recall the procedure in the special case of (5) when 
every r; = 1 (i.e., ordinary rational interpolation) and each a; is found from 


) — Piles) _ Gi0pea(as) + (ei — Tea)pi-2( zi) 
(6) Kes) qi(x4) 4i,0Qi-1(24) + (% — Lis) Qi-2(i) ‘ 





where p;(x) and g;(z) denote the numerator and denominator of the ith convergent 
[2]. (We may differ here from some standard notation by referring to a; in (5) as 
the first, instead of the zero-th convergent. ) 

The same idea as in (6) can be used for the stepwise determination of the entire 

















block of r; coefficients ajo, @i1, --* , @,r,-. simultaneously. Assuming that from 
knowledge of ai0, 41,1, ***, @i-1,7;_,1 we have found in (5) the convergents 
Poa(%)/Qea(2), De-2(t)/qe-o(z), Where s — 1 = Dost r;. Then 
Ri(x)paa(z) + (x — t1)po(z) 
7 N D(z) = ' 
(7) (2)/D(=) = “Pa)q.-a(a) + @ — 2eaqe-al@) 
where 
Rx) = aio +> Sel gp iF ses 
| Gi | Qi,2 

(8) 

42 —%| zr—a| , t— in| pha tel 

| irs i+1,0 | Gi41,1 | Ga.rn—l 
We may neglect that part of R;(«) which proceeds from 7 ~ x: | + .-- in determin- 
++1,0 

ing Gio, @i,1, -** , @,,;-1 from the r; osculatory interpolating conditions on (7) at 


x = 2x;. Denoting that truncated R,(x) by R;*(x) = S.(x)/T(x), we apply the 
general lemma at x = 7; , finding first the coefficients of the polynomials S;(x) and 
T (x) to satisfy 
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{Si(x)pea(xz) + (2 — 2-1) Tx) pro(x)}™ lone, 
(9) = [f(x){ Six) qea(z) + (x — 2a) T(x) qe-2(z)}) | ene, 5 
m=0,1,---,r—1. 


The next operation is to retrieve the coefficients ajo, --- , @:,7,-1 from those of S,(x) 
and 7'(x). Finally, the convergents to (5), p.(x), gs(%), Deri(Z), Ge4i(%), --- are 
found by the usual recurrence scheme 


fort = 0 


(x — Xi1) ’ 
for ¢ = 1(1)r; — 1, 


(10) Port(Z) = i,t Dett—1(2) + (x — 2;) 


Ps+t—2( x ), 


(10’) Gere( 2) = i,t Jort—a( 2) + 4 nar Qe+t-2(2), pn = a ~~ 
The cycle of (9), retrieval of coefficients a;,_ and (10), (10’) is repeated, but with 7 
replaced by i + 1 and s replaced by s + r;, etc. For a succession of smaller values 
of r;, say 7; = 2, riz, = 3, --- it may be convenient to apply (9) at several points 
at once. 

This procedure, which involves the solution of a number of separate linear sys- 
tems of orders r;(or 7; + rig: + --- + ri4; when fairly small), is an intermediate 
one between that of solving for all >>2., r; coefficients of N(x)/D(x) at one time, 
where the linear system may be inconveniently large, and that of solving for each 
separate coefficient a;,, in (5), from an equation like the following, when m > 0, 


ae Co oo im Dr(X) + (x en speis)}" 
(iL) f(x) jean) re gern 





i-—1 
’ f= > T; + m, 


rz; j=1 





Scuepute I: S,(x)/T (zx) 








ri S,(x) T (x) 

1 | a 4 y 
2 | ayao + (x — 2;) ay 

3 | G2Qyd0 + (a2 + ao) (& — 24) 2a, + (x — 24) 


4 | Q3Q2010 + (asQ2 + a3Qo + a0) (x — %4) | @3Q2a, + (az + a1) (z — 2,) 
+ (x — 2;)? 


5 | GgQgQ201Ao + (A432 + ayAzAo 43020, + (aya3 + aya; 

+ a4Q,Q9 + 20109)(x — x,) + a2a;)(x — x4) + (x — 2)? 
+ (as + a2 + ao)(x — 2,)? 

6 | GgQgQgG2Q1A9 + (AsQgAghe + Ag Q43Qo | AsQgQ30201 + (sQ4Q3 + A504: 

+ 504010) + A5201Ao + A3020,M)(x | + ds@201 + a34201)(x — x;) 

— 2) + (G54 + age + asQo + 302 + (ds + a3 + a)(x — 2;)? 

+ a3o + a)(x — x;)? + (x — 2;)3 














EE 
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(remembering now that we must retain the denominator in the right niember of 
(11) because it is for a single m) and finding N(x) and D(z) from (10), (10’). 
Since there would be >>?.; r; such equations (11), that is likely to involve more 
work than the present scheme. 

In using (9) at just a single point z;, S;(z) and T,(z) are kept in terms of 
(x — 2)’, whereas p(x), qe-1(1), Pe-2(z) and Ge-2(x) are in terms of x’. After 
retrieving the coefficients a;,,,, it is natural to obtain p,4:(x), ge4:(z) from (10), 
(10’) in terms of x’, then S;;(z) and 74:(z) are found as polynomials in x — 2j4;, 
etc. 

For the convenience of the user we list in ScuepuLE I the explicit expressions 
for S;(x) and T(x) as polynomials in (x — z;,), for r; = 1(1)6, and the sequence 
of operations for retrieving the coefficients a;,.,. For the sake of brevity we write, 
from now On, dm fOr Ai,m - 

Denoting the coefficients of (x — 2z;)’ in S,(x) and T(x) by a; and 8; respec- 
tively, the coefficients a,, are retrieved as follows: For r; = 1 and 2, obvious. For 
T; = 3, do = a/Bo, Ge = a — Ao, 0, = Bo/ae. For r; = 4,09 = ao/Bo , a1 = Bo/(an — 
of), a3 = By — Gy , G2 = (a — Oef)/az. Forr; = 5,09 = a0/Bo , a1 = Bo/ (a1 — ahs), 
a2 = (a, — aofs)/{Bi — ay(a2 — do)}, @s = ae — Go — G2, a; = [B, — a( — a)} 
/ay. For ri = 6, do = ao/Bo, a: = Bo/(a1 — afi), a2 = (a1 — as)/{fi — 
a;(a2 — aoS2)}, a3 = 18: — ai(a2 — aq2)}/{a2 — (ao + ae)(B2 — ai) — ayo}, as = 
Be — G3 — G, 4 = {ar — (do + a2)(B2 — a1) — ayao}/as. 

The retrieval process assumes that every a,, ~ 0 for m ~ 0. But, we may have 
a = 0, and every retrieval step still goes through. 

It seems likely from the above (verification left as an exercise for the reader) 
that for any r; there is a general procedure for retrieving ao , a: , --* , @r,—3 5 @r;-1, 
a,,-2, m that order, from the Euler-Minding formula [3], according to which 
S,(z) and 7;(z) have these explicit expressions: 

0,r;—2 


Sx) = aoq, --- Oy ;-1 (1 + Zz (z — Li) / 05; Aj, 


0,r;—3 


(12) + > ® (z — 24)" /05 0541 0n4 Ar+2 
j<k 


3 


0,r;—4 
+ > pape (x — 2¢)°/aj 0541 0n41 04420142043 + ), 
<k< 
and 
1ri-—2 
T(x) = GiG@e °°- an. (1 + > (x — Xi) /Q; A541 
(13) ‘ 


1ry-—3 
+ 0D (x — 25)*/aj 0541 0e4Oey2 + +), 
j<k 


General Dynamics/Astronautics 
San Diego, California 


1. L. M. Mrtnz-Tuomson, The Calculus of Finite Differences, Macmillan, London, 1933, 
Chapter V, p. 104-123. 

2. P. I. Ricnarps, Manual of Mathematical Physics, Pergamon, London and New York, 
1959, p. 257. 
3. O. Perron, Die Lehre von den Kettenbriichen, B. G. Teubner, Stuttgart, 1954, Vol. I, 
p. 5-7. 
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On Solving y’ = fy + g with a Boundary 
Condition at Infinity 


By Charlotte Froese 
Consider the class of differential equations 
(1) y” = f(z)y + g(x) 
for which 


(i) g(t) -0 as tr @ 
(ii) f(z) > 0 and f(t4)—c¢ as tr—> ~, 
together with the boundary conditions, 


(ili) y(a) = yo and y() = 0. 


Suppose further that a numerical solution is required over the range (a, b) where 
b is such that | y’(z) | < +r, for all z = b. 

For large values of z, one solution of the equation tends to infinity whereas the 
other tends to zero; the latter is the desired solution. Because of round-off and trun- 
cation errors which are inherent in a numerical procedure, any outward integration 
will introduce some of the former solution into the calculation, so that the numeri- 
cal solution will not tend to zero but eventually increase exponentially. This diffi- 
culty may be avoided by determining the solution at some large value of z from an 
asymptotic expansion, for example, and integrating inwards. It is avoided more 
simply by the procedure to be described, particularly when the inward integration 


cannot be started readily and the range (a, b) over which the solution extends is 
not known in advance. 


Let 
Xi = Xo + th, X= a, 


and define y; to be the computed approximation to y(x;); yo is known and 4; , ye - -- 
are to be determined. 
It is well known that 


2 
(2) ay; =} (ue + 5 v’) + 0(h°). 


Substituting for y” from (1) and neglecting terms which are 0(h°), we may rewrite 
(2) as a system of equations 


—bi a 

Ce oe eae 
a - 

| Prenton \o cs 
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where 
2 
a;= 1 a7 
i> 2+ ORY 
and 


2 
= - (gisr + 10g: + gi), 


or, in matrix form, AY = C. The matrix A has a tridiagonal form so that the system 
of equations may readily be solved by the Gauss elimination method. An algorithm 
in this case is given by equations (3) and (4): 


d, = —h d; = —b; — (a;4a;/d;4) 


Za = C1 — Ao 25 = C5 — (A:-42i1/di4) 


(3) 


and 


(4) Ys = (25 — GigsYizs)/di 


Unlike the standard case, the order N of the matrix is not known in advance but 
depends on the solution. For f(z) > 0 and h sufficiently small, it can be shown that 
d; — —1 + O(h) as i increases and obviously a; = 1 + 0(h’). Thus, for i suffi- 
ciently large and h sufficiently small equation (4) becomes 


Yin — Yi ™ 2; 


or 

a™~ hy. 
Therefore, if we solve equations (3) fori = 2,3, --- , N successively until | zy | < hr 
followed by (4) withi = N, N — 1, --- , 1 we will have obtained a solution over 


the desired range. The solution of (4) when i = N requires an approximation; the 
simplest is 


(i) yw = 0. 
A more accurate assumption is 
(ii) ywar = Cyn, 
where C is determined from an asymptotic expansion. Then 
Yn = 2n/(dw + Gwu:C). 


The above boundary conditions are frequently associated with eigenvalue prob- 
lems. An example is the Hartree-Fock equation for an electron, 


~_f .2¥@), +D\,_ 
P -(. 220) 4 WED) p - feor. 


Here ¢ is an eigenvalue and a solution of the problem may not exist for an arbitrary 
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choice of «. The usual procedure is to divide the range of integration into two parts, 
integrate outwards for a solution satisfying one boundary condition, integrate in- 
wards for a solution satisfying the other boundary condition, match the solutions at 
an intermediate point and adjust « so that the derivatives also agree [1], [2]. The 
inward integration may be avoided with the procedure described earlier. A con- 
venient way of dividing the range is according to the sign of f(r). For some r, 
f(r) < 0so that condition (ii) is not satisfied: the procedure described here is not 
always numerically stable when f(r) < 0 [3]; in fact, for some values of i, | d; | < 1. 
Of a series of standard methods, the Numerov method, 


2 
Ynt+i — ((2 + s itn) Yn rt (1 _. ¥ fu+) Yn-1 


+ - (gn+1 + 10gn + m))/(1 = fons) ’ 


was found to be most accurate in this case, for a given number of evaluations of f 
for the outward integration. The procedure used successfully was to integrate out- 
wards according to (5) until f(r) > 0, then, with the last value computed as a boun- 
dary condition, to solve for the “tail” of the wave function by the method described 
here. The energy adjustment will be the same as before. 

Department of Mathematics 


University of British Columbia 
Vancouver, B. C. 


(5) 


1. J. W. Cootey, ‘‘An improved eigenvalue corrector formula for solving the Schrédinger 
equation for central fields,’’ Math. Comp., v. 15, 1961, p. 363-374 
2 


. A. 8. Dovetas, “On the Sturm-Liouville equation with two-point boundary condi- 
tions,” Proc. Cambridge Philos. Soc., v. 52, 1956, p. 


3. E. C. Riptey, “A numerical method of cadens second-order linear differential equa- 
tions with two-point boundary conditions,’’ Proc. Cambridge Philos. Soc., v. 53, 1957, p. 442-447. 


On the Inversion of Sparse Matrices 


By A. L. Dulmage and N. S. Mendelsohn 


1. Introduction. There are a number of problems in applied mathematics involv- 
ing many equations in many unknowns, but for which each equation involves only 
a small fraction of the unknowns. If such problems are linear or are approximated 
by linearization, one is involved with a matrix, a large proportion of whose entries 
are zero. To invert such a matrix A it is sometimes advantageous to permute the 
rows and columns of A yielding PAQ where P and Q are permutation matrices. If 


| Ay 0 | 
| As | 
PAQ = | —— 
| a 
wy A, | 

where A; , A>, --- A, are square matrices, the problem of inverting PAQ is reduced 


Received September 18, 1961. This research was supported by the United States Air Force 
Office of Scientific Research. 
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to inverting the blocks A, followed by matrix multiplication and addition. On per- 
muting the rows and columns of (PAQ)~™ one obtains the matrix A™. 

F. Harary [4] gave a method based on the connectivity theory of directed 
graphs. The blocks A; turn out to be the matrix representatives of the strong com- 
ponents of an associated directed graph. Harary’s method requires that Q = P™’,a 
restriction which is quite unnecessary for matrix inversion. As a result, many 
matrices which reduce under independent permutations of rows and columns will 
not reduce if one insists that Q = P™*. To remove this latter restriction, the authors 
have replaced a directed graph by a bipartite graph. The strong components of a 
directed graph become the irreducible components of a bipartite graph (see [1] 
for the definition of an irreducible component of a bipartite graph, [2] for the con- 


nection between strong components of directed graphs and irreducible components 
of bipartite graphs). 


2. The Method of Reduction. Let A be a square matrix of order n with entries 
a;;. Associated with A is a bipartite graph K, with two sets of vertices S = s,, 
S,-::,8, andT = 4,4, ---,t,.A pair (s;, t;) is an edge of K, if and only if 
ai; ~ 0 (one obtains Harary’s directed graph if one identifies the vertices s; and 
t:). 

Suppose we can find matrices P and Q such that 


| Ay 0 
Zz 
PAQ =| 
* A, | 
and such thai A; , Az, --- , A, cannot be further reduced by permutations of their 
rows and columns. Then the graphs corresponding to A; , Az, --- , Ar are the ir- 


reducible components of K, . 

In [3], the authors give an algorithm for obtaining the irreducible components of 
a bipartite graph. This algorithm is easily programmable for machine computation 
and yields the permutations P and Q. In the case where A is a non-singular matrix 
the graph K, has no tails (see [3] for definition of a tail) and the algorithm de- 
scribed in [3] can be considerably simplified, since the steps needed to isolate and 
identify the tails can be omitted. An alternative procedure is the following. First, 
locate a nonzero term in the expansion of the determinant of A. This can readily be 
done using the algorithm of Marshall Hall [5] or Fulkerson and Ford [6]. Next, 
permute the rows of A until the entries of this term occupy the main diagonal. Call 
this new matrix A*. Finally, apply the method of F. Harary given in [4] to A*. 

We append an example: 


lett A= 








—_wworods 
onorono 
corocooonw 
ROoOwWNSOOS 
onwooceo 
conocer 
cocoooro 
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Note that the corresponding directed graph has the Hamiltonian circuit 1 — 6 — 
5—3—4-—2-—7 — 1 and so is strongly connected. Hence, for no permutation 
P does PAP™ reduce. Using the algorithm described in [3] one obtains the permuta- 
tions P = (1725643) and Q = (1)(236754). Applying P and Q to the 
rows and columns of A, one obtains: 











1 210 0 O10 0 
14,0 0 O10 0 

0 314 0 O10 O 

PAQ =|3 O0j2 2 O}0 O 

0 012 0 4;0 0 

3010 3 O11 2 

0 010 0 015 1 

The authors are indebted to F. Harary for a pre-publication copy of his paper [4]. 


The University of Manitoba 
Winnipeg, Canada 


1. A.L. DULMAGE & N.S. MEenpEtsonn, “‘A structure theory of bipartite graphs of finite 
exterior dimension,”’ Trans. Roy. Soc. Canada, Sect. III 53, 1959, p. 1-13. 


2. D.M. JoHNSON, A. L. Dutmageg, &N. 8. MENDELSOBN, “Connectivity and reducibility 
of graphs’ Canad. J. Math., 14, 1962, p. 529-539. 
3 


Dutmace & N. 8. MENDELSOEN, “‘Two algorithms for bipartite graphs,’’ to be 
eS in Soc. Indust. Appl. Math. 


Harary, “A graph theoretic approach to matrix inversion by partitioning,’’ to be 
4 & in Numer. Math. 


M. Haut, ‘‘An algorithm for distinct representatives,’’ Amer. Math. Monthly, 1956, 
p. 716-71. 


. L. R. Forp & D. R. Futxerson, “A simple algorithm for finding maximal network 
fous ‘and an application to the Hitcheock problem,” Canad. J. Math. 9, 1957, p. 210-218. Also 
in Management Sci., October 1956, p. 24-32. 


Missing Data Correlation Computations 


By R. L. Jennrich 


In correlation analysis or in any multivariate analysis based on the computation 
of a correlation or covariance matrix, the applied statistician often runs into the 
problem of missing data. To avoid complication in computing the correlation 
matrix, a complete observation vector is often discarded when only one or more of 
its components are missing. If a correlation matrix is computed by means of a 
standard electronic computer program, this procedure is often necessary. A large 
percentage of data may be thrown away when only a small percentage is missing. 
This note describes a modification in the standard computing scheme which elimi- 
nates this waste of data. 

Let 2n1 , Zn2, *** » Lap denote the p components of the nth observation vector, 
n = 1,2, --- ,N.Itis customary to add an n + Ist component to this vector which 
is identically equal to one. That is 2,54: = 1. The cross product matrix 


N 
Qj = Ditates §=1,-+-,p4+1; j=1,---,p+1 
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is customarily computed and the correlation matrix is computed from this matrix 
by using the fact that 


N 
Orit.ets = N, Giees = Ltt; t= 1,---,p. 
n=l 


In addition to adding a component which is identically one to each observation 
vector, let us form a new vector Cn: , Cn2, *** , Cap Where C,; is zero if the ith com- 
ponent of the observation vector is missing and one otherwise. Letting each element 
of missing data have value zero, we form the cross product matrices 


N 

8ij = Do tains i_j=i,---,pti 
N 

mis = Do nitns ,j=1,---,p. 


The means m;, covariances v;;, and correlations r;; are computed from these 
matrices by the formulas 


1 


mM; = — 8i,ptt 
Ni 
oo. Fe aay 
Ny 
Vij 
Ti = 


It should be noted that the statistical properties of these estimates will differ 
slightly from those computed without missing data. A discussion of some of these 
properties is given by 8S. 8. Wilks [1]. 

A FORTRAN program for the computations described in this note is in use at 
the University of Wisconsin. A write-up and program deck can be obtained by 
writing to the author. 


The University of Wisconsin 
Madison 6, Wisconsin 


1. S. S. Wixxs, ‘“Moments and distributions of estimates of population parameters from 
fragmentary samples,”’ Ann. Math. Stat., v. 3, 1932, p. 163. 


Polynomial Approximations to I(x), I,(«) and 
Related Functions 


By F. D. Burgoyne 


Hitchcock [1] gives polynomial approximations to some Bessel functions of order 
zero and one and to some related functions. Notable omissions from his list are any 
approximations to Io(z) or I,(z). The following approximations may serve to fill 
this gap. 

If we write I,(z) = (2er)~*e*F,(x), then with the maximum error stated in 
brackets in each case, and provided 0 < ¢ < 1, 
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I(4t) = 0.99999 99985 + 4.00000 01935 ¢ + 3.99999 59541 ¢ 
+ 1.77780 99690 ¢° + 0.44431 89384 ¢° + 0.07137 58187 ¢”° 
+ 0.00759 42968 ¢” + 0.00082 67816 ¢* (17 x 107”), 

£"I,(4t) = 1.99999 99997 + 4.00000 00421 ¢ + 2.66666 57853 ¢ 

+ 0.88889 59049 ¢° + 0.17775 04042 ¢® + 0.02376 15011 ¢° 
+ 0.00219 03549 #” + 0.00020 11611 7“ (4 x 10°”), 

(2x) "*F,(4/t) = 0.39894 22809 + 0.01246 67783 t + 0.00176 23668 ? 
+ 0.00026 22220 ¢ + 0.00225 85672 ¢# — 0.01283 14822 ¢ 
+ 0.04958 11198 ¢° — 0.12099 40805 ¢’ + 0.18954 76618 ¢° 
— 0.18677 83276 # + 0.11133 15511 ¢° — 0.03666 94167 
+ 0.00512 46015 ¢” (7 x 10°”), 

(2) *?F,(4/t) = 0.39894 22799 — 0.03740 06642 ¢ — 0.00293 14981 ¢ 
— 0.00043 77220 # — 0.00237 87859 ¢ + 0.01319 50213 ¢ 
— 0.05078 72951 t° + 0.12301 43060 ¢’ — 0.19083 32956 ¢° 
+ 0.18552 23758 # — 0.10862 98349 ¢° + 0.03497 54315 7? 
— 0.00474 86397 #” (8 x 10°”). 


The first two approximations were obtained by the economization method of 
Lanczos [2], which is used by Hitchcock. As he notes, this method is inapplicable 
for the last two approximations, and these were obtained by collocation at the 
zeros of TT;(x) = cos {13 cos (2x — 1)}. 


Battersea College of Technology 
London, 8.W. 11. 


1. A. J. M. Hircucock, ‘“‘Polynomial approximations to Bessel functions of order zero and 


one and to related functions,” MTAC, v. 11, 1957, p. 86-88. 
2. C. Lanczos, Applied Analysis, Prentice Hall, Inc., New Jersey, 1956. 


A Note on the Curve Fitting of Discrete Data 
by Economization 


By F. D. Burgoyne 


Suppose that we are given a set of points (z;, y;) 0 S i S n and we desire to 
find the polynomial p(x) of given degree m(<n) such that max; |! y; — p(z,) | is 
a minimum. It is well known that this may be performed in good approximation 
by using the method of least squares to find the polynomial g(x) of degree m such 
that >>: {ys — q(z;)}* is a minimum, and then taking p(z) = q(x) + c, where ¢ is 
constant given by 


2c = min; {y; — g(x:)} + max; {y; — ¢(z,)}. 
Received July 13, 1961. 
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While several methods (mainly of an iterative nature) exist for finding a better 
approximation to p(x), e.g., that of Hastings [1], none seem quite so simple as the 
above. However, the author has for some time also used the following procedure 
which is both simple and straightforward and does not involve iteration; it is in 
fact basically the analogue for discrete data of the usual economization process. 
First, the polynomial of degree n is found which passes exactly through the points 
(xi, yi), then this polynomial is economized to a polynomial r(z) of degree m over 
the range min; 7; S x S max; z;. We then take p(z) = r(x) + d, wheredisa 
constant given by 


2d = min; {y; — r(zi)} + max; {ys — r(x}. 
When n is not too large and the z; are equally spaced r(z) may be found by 


hand as follows. First the z; are transformed so that z; = i/n. Now by Newton’s 
interpolation formula the polynomial of degree n which passes exactly through the 


points (2;, ys) is Dx 7 A‘z), and the polynomial raim(2z), which (7) be- 


comes when economized to degree m over the range 0 < xz S 1, may be found 
from a previously prepared table. Then r(z) = do: Taim( 2) A’zy . If desired, a 
different interpolation formula may be used. 

In the vast majority of examples tried by the author it was found that 


#max; | ys — 9(2:) —ec|< max; |y; — r(z;) — d| S max; | y; — g(x) —c|. 


In each of the few cases in which max;|y;— r(z;) — d| was greater than 
max; | yi — g(z:) — c| it was found that the points (z;, y;) gave an inadequate 
picture of the polynomial of degree n which passed exactly through them; this 
rarely occurred when the xz; were equally spaced. 

For some of the examples p(x) also was found by using Hastings’ method [1]. 
In each of these cases max; | y; — r(xi) — d| was less than $max; | y; — p(z,) | . 

As an illustration consider the following example. We are given the points 
(0, 1), (0.25, 1), (0.5, 15), (0.75, 79), (1, 253), and we desire the quadratic poly- 
nomial which best approximates them. In this case p(x) can be found exactly, and 
we obtain the following results: 


p(x) = 17 —2442 +4642", maximum error 16; 
r(x) +d = 15.5 —2282 +4482", maximum error 17.5. 


The maximum error associated with q(x) + c is 19.2. 

It may be added that on Pegasus both the least squares and the economization 
fitting processes take roughly the same time and use roughly the same number of 
storage locations. 


Battersea College of Technology, 
London, 8.W.11. 


1. C. Hastines, Approximations for Digital Computers, Princeton University Press, 1955. 
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55[F].—Kenneta I. Appet & J. Barxiey Rosser, Table for Estimating Functions 
of Primes, Communications Research Division Technical Report No. 4, Institute 
for Defense Analyses, Princeton, N. J., 1961, xxxii + 125 p., 22 cm. 

This interesting table is quite complicated in its design, and partly so in its sub- 
ject matter. For the sake of brevity a complete description will not be given here. 
The functions of the primes referred to in the title are primarily r(x), the number 
of primes <2, the three sums over the primes p: 


1 ] 
> log P; > on pm =ae, 
pst pszP psz Pp 


and the product 


mo 
psz Pp — 1 
The range of z is from 2 to 10°. 
The approach used here may be illustrated with the first sum: 


6(z) = > log p. 
Pst 


It is well-known that 0(2) ~ 2, (this is equivalent to the Prime Number Theorem), 
and if one defines the coefficient TH(x) by 


(2) = x — TH(z)-vz, 


it is found empirically that TH(x) has a rather small variation over the range x 
for which it has been computed. For example, between zt = 85,881,353 and z = 
87,679,913, TH(z) has a maximum of 1.337 and a minimum of 1.015. Knowing 
these extremes, one could therefore obtain bounds for, say, @(86,692,297), as 


follows: 
86,679,848 < 6 <= 86,682,847. 
The exact value for this argument is @ = 86,681,759.3. 
Similarly, if one knows bounds on 


PI(x) = [li(x) —x(zx)] log x-2”, 


SR(x) = x" logz [> ; — log log z — 0.261497 | : 


Pst 


and 


SL(z) = 2” [> wee — logz + 1.332582 |, 
psz 


one can obtain good estimates for x(x) and for the second and third sums above. 
Further, the authors also give more complicated, third-order correction coeffi- 
cients which enable one to obtain even more accurate estimates. 
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The corresponding coefficient for the product: 


PR(z) = ole II My cs log s| ; 
pazp—1 
is not listed as such, since the authors, by study of their preliminary computations, 
discovered the approximate formula: 


a (SR(x))* log x 
PR(z) r= SR(2) + oomings — (1 + loge) 
Professor Lowell Schoenfeld later obtained a rigorous justification. His long proof 
is given in full in the Introduction. 

The five primary functions are list.d, first, for each of the first 64 primes; x(x) 
exactly, 0(z) to 5D, and the other three functions to 10D. The primes from 313 to 
99,999,989 are divided into 173 rather irregularly spaced intervals, and in each 
interval the five functions are listed at every z for which one of the auxiliary func- 
tions, TH (x), PI(xz), etc., takes on a maximum or a minimum value. These extrema 
of TH(z), ete., are also given. 

In addition, the largest gap between successive primes in each interval is listed. 
The largest gap up to 10° occurs between the primes 47,326,693 and 47,326,913. 

The orientation here is that of estimating the primary functions in terms of the 
bounded coefficients. But from a theoretical point of view the opposite orientation 
is probably one of greater interest. One would like to know the order of the true 
bounds upon these coefficients. The z"” that enters into all of their defining equa- 
tions is related to the Riemann Hypothesis, and, as is well-known, the state of the 
theory here leaves much to be desired. It has been suggested, in an off-hand manner, 
in MTAC, v. 13, 1959, p. 282, that PJ(z) has a mean value equal to 1. The range 
of its values given here, for 313 < x < 10°, is 


0.526 = PI(x) S 2.742. 





Finally, a word concerning the table per se. Since the subject matter is so funda- 
mental, an improved and more elegant edition is probably called for. While the table, 
as it stands, is quite workable, a less erratic selection of intervals and a somewhat 
clearer format would be desirable. 


D.S. 


56[G].—G. E. Suitov, An Introduction to the Theory of Linear Spaces, Prentice-Hall, 
Inc. New Jersey, 1961, ix + 310 p., 23 cm. Price $10.00. 


This book is the first in Prentice-Hall’s series of translations from the Russian. 
A bibliography has been added by the translator, R. A. Silverman. 

The contents include the usual topics in linear algebra such as determinants, 
linear spaces, systems of linear equations, coordinate transformations, invariant 
subspaces and eigenvalues of linear transformations, and quadratic forms. The 
degree of abstraction is shown by the fact that sections on ideals and tensors are 
included, but marked with asterisks to indicate that they may be omitted if desired. 
The final chapter deals with infinite-dimensional Euclidean spaces. 

The book is not concerned with computing methods directly, so that its value 








502 REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


to one interested in numerical work lies in its development of needed topics in the 
theory of matrices. Its lucid treatment of the topics covered makes it a fine addition 
to the literature. The inclusion of suitable problems also makes it useful for class- 
room use. 

Ausert A. Grau 


Oak Ridge National Laboratory 
Oak Ridge, Tennesee 


57[J].—_L. B. W. Jotuxy, Summation of Series, Dover Publications, Inc., New 
York, 1961, xii + 251 p., 20 cm. Price $2.00 (Paperbound). 


This is a “revised and enlarged version of the work first published by Chapman 
& Hall, Ltd. in 1925”. The 700-odd series in the former edition have now been 
increased in number to 1146. For most of these series a specific reference is listed. 
While there is much of use and interest here, there are also, in the opinion of the 
undersigned, numerous defects. 

The notation used is often disturbingly original. Thus: 





(71) 1—-4F4+4-444--- © = logh2 
(97) \eerpee yes, ene 
(361) oe =5, 

(373) > aay = # 


(1133) Sum of Power Series 
1 1 1 
oe ak os Se 2 


There are misprints and resulting obscurities: 


B.. Ph, 1 oleh 
(94) Be =1+stst Reo ee 
(335) if >> = = {(s) 2,3,5 --- p—are prime numbers in order 
1 


ts + 38+ 38+ --- © = ¢(s)(1 — 2°) 


It is not clear what is being “summed” in: 


(68) 2+5+13+ 35 + --- n terms = }(3" — 1) +2"-1 
(793) @ — 26° + --- © = logh (1 + Osin 0) 

a 7 24 = 
(794) — 3 oe! — °° @ = logh d cot 6 


since the continuations of the series on the left are not at all obvious. 
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The vigorously divergent 


(976) 1 (1 = 214.31 — +++ @) = SAO 
is given without any tiresome commentary concerning convergence. 

While a good handbook of infinite series is something much to be desired, it is 
doubtful whether the present book fully meets this need. 


D. 8. 


58[K].—B. M. Bennett & P. Hsu, Significance Tests in a 2 X 2 Contingency Table: 
Additional Extension of Finney-Latscha Tables, March 1962. Deposited in 
UMT File. [See Review 9, Math. Comp, v. 15, 1961, p. 88-89; Review 20, ibid., 
v. 16, 1962, p. 252-253.] 


The authors present in these manuscript tables still another addition to the tables 
of Finney, first extended by Latscha. This latest extension covers the range A = 
31(1)45, giving the exact test probabilities to four decimal places, as previously, 
and the significant values of b at the four levels presented in the Finney-Latscha 
tables and retained in the previous extensions thereto by the present authors. 


J. W. W. 


59[V].—Puiir D. THompson, Numerical Weather Analysis and Prediction, June 19, 
1961. The Macmillan Co., New York, xiv + 170 p., 24 cm. Price $6.50. 


In a rapidly developing field it never seems quite appropriate to freeze the state 
of knowledge in the form of a book. However, enough has evolved since the begin- 
nings of numerical weather prediction to warrant a knowledgeable appraisal of the 
course of its development. Not only should such a book have a didactic objective 
but one would hope that the perspective be equally useful as a reference for active 
workers in the field. This would require the text to assess the road of experience 
well enough to define the problems and to indicate the avenues which are likely to 
yield a fruitful expansion of knowledge. Colonel Thompson’s book represents a first 
such attempt. The fact that he has contributed materially to the evolution he sets 
out to document, taken together with his characteristically smooth expository style, 
amply qualify him for the task. 

In Chapter 1, after a brief discussion of the inherent difficulties of observing and 
forecasting the atmosphere’s evolutions, one finds a description of instrumental, 
aerological, and analysis techniques, and of the atmosphere’s kinematical character- 
istics. The author then indicates the role of hydrodynamic laboratory models as a 
research tool and gives an historical development of numerical weather prediction: 
the antecedents heralding the Norwegian school and the contributions of Richard- 
son, Rossby and the Princeton group. 

Chapters 2 and 3 are given over to a summary of the hydrothermodynamic 
equations of meteorology, first in height coordinates and then in terms of pressure 
and of potential temperature. The transformation and interpretation of upper and 
lower boundary conditions are not given. Methods of central differencing and practi- 
cal problems of numerical weather prediction, especially those engendered by the 
quasi-non-divergent character of the atmosphere, are then discussed. Chapters 4 
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and 5 first go into the properties of pure sound, gravitational, and Rossby waves by 
means of a perturbation analysis of the linear hydrodynamic equations, followed 
by an analysis of the corresponding finite difference equations with some discussion 
of their computational stability properties. As an introduction to the use of filtering 
approximations, the author analyzes in Chapter 6 the properties of mixed wave 
type solutions and their interaction in linear systems. 

In Chapter 7 the author discusses physical systems in which the total kinetic 
energy is preserved by means of the equivalent barotropic approximation, giving 
some emphasis to the numerical solution of the corresponding non-linear equations. 
Chapter 8 lightly covers the question of mapping, examines the finite Fourier series 
method of Charney, Fjortoft, and von Neumann, and then goes on to discuss the 
application of relaxation techniques to elliptic difference equations. 

Systems in which potential-kinetic energy conversions are admissible and the 
necessary vertical differencing structure are considered within a geostrophic frame- 
work in Chapter 9. Here only superficial mention is made of the approximations 
energetically consistent with the geostrophic constraint. Furthermore, the question 
of the consistent lateral boundary conditions for the vertical velocity and thermo- 
dynamic equations is not covered at; all. In Chapter 10 Dr. Thompson analyzes the 
process of baroclinic instability, and the concomitant energy conversions and pole- 
ward heat transfer, and their role in the index cycle. However he does not take the 
opportunity to indicate how the barotropic kinetic energy of vortices is transformed 
to maintain the westerlies—the final link in the energy cycle. 

The subject of Chapter 11 is the balance approximation as a filtering device. It 
could have been useful here to show its connection with the primitive barotropic 
equations. Also absent is mention of attempts to adapt balance techniques to baro- 
clinic models. In Chapter 12 he takes up the question of establishing initial condi- 
tions for the primitive equations. He apparently holds the opinion that the balanc- 
ing constraint must be applied periodically; however, this has not been borne out 
by experience. He considers the question of formulating boundary conditions for 
open boundaries as an unsolved problem. 

Getting down to practical problems, the author discusses in some detail in 
Chapter 13 the question of operational utilization. In particular, he describes the 
organization of the Joint Numerical Weather Prediction Unit, methods of data 
processing and objective analysis. After drawing comparisons with routine conven- 
tional forecasting methods and skills, he then boldly makes an attempt to calculate 
the economic worth of a forecaster. Finally, he examines the impact of mechaniza- 
tion on the data processing chain. 

Chapter 14 is given over to the question of unsolved problems. Here he points 
out the systematic errors in the zonal angular momentum which result from ignor- 
ing surface stresses. He then goes on to emphasize the importance of sources of 
kinetic energy through baroclinic instability but makes no mention of how the 
available zonal potential energy is maintained through radiative processes, which 
are probably of some significance within 48 hours. He makes some mention of how 
the quasi-stationary long waves may be excited geographically and of the uncer- 
tainties of the initial state on a forecast. Finally, he says something about the 
effects of truncation and round-off error. Chapter 15, entitled ‘““The Outlook for the 
Future,” is essentially a recapitulation of the 14 preceding chapters. 





a» lf es Go te A ce oe oe ae me lL Ue 6 
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One’s impression upon reading the book is that it is more an essay than a text or 
reference. A good many of its 170 pages are given over to the discussion of peripheral 
subjects which are treated much more exhaustively elsewhere and could more 
appropriately have been referred to. This reviewer would rather that the author 
had taken this space and perhaps more for greater thoroughness in discussing prob- 
lems intimately germane to the development of dynamical prediction by numerical 
methods. A more complete discussion of computational instability of the various 
types that have already been encountered would have been extremely useful. A 
comprehensive account of mapping techniques for finite differences would be useful 
if found in one place. Very little attention is given to Green’s function techniques 
and Fourier space, to the process of barotropic stability, to Lagrangian methods, 
to “staggered’’ finite difference methods yielding non-redundant solutions (such as 
that of Eliassen). The powerful methods and useful results of scale analysis are 
prominent in their omission. One would hope that in the absence of thoroughness 
Dr. Thompson would have given an exhaustive bibliography; however, his refer- 
ences are sometimes vague if not scanty (Phillips’ significant contributions are 
virtually ignored) and those which are included often are given a superficial critique. 
Being the first of its kind, this book does fill a gap. However, this reviewer feels it to 
fall short of the needs, if not of the author’s objectives. 

JOSEPH SMAGORINSKY 


U. S. Weather Bureau 
Washington 25, D. C. 


60{X].—_W. F. Frerpercer, Editor-in-chief, International Dictionary of Applied 
Mathematics, D. Van Nostrand Co., Inc., Princeton, N. J., 1960, 1173 p., 
26 cm. Price $25.00. 


At the outset it is perhaps appropriate to say a word concerning the title of this 
useful reference book. One might think that a book so named would confine itself 
to descriptions of those branches of mathematics which are applied to physics, 
engineering, etc., that is, to numerical analysis, vector analysis, statistics, etc. 
Instead a large number, or even most, of the entries here, e.g., Binary Stars, Polymer, 
Isotopes, Pfund Series, etc., are descriptions of those phenomena to which such 
mathematics may be applied. Of the 32 fields covered in this volume only 6 are 
applied mathematics in the strict sense, while Acoustical Engineering, Acoustics, 
Aerodynamics and Hydrodynamics, Astronomy, Atomic Structure, Automatic 
Control, Chemistry, Elasticity, Electromagnetic Theory, and 17 other fields are, 
rather, physical sciences to which mathematics is applied. 

Each of the 32 fields had one or more authorities as a contributing editor. For 
example, that for Numerical Analysis was A. 8. Householder. The 8000-plus entries 
are all listed alphabetically and not by field. 

The many entries differ widely in their length and character. Those on modern 
physics, e.g., Relativity, Quantum Mechanics, S-Matrix, etc. are often fairly long 
and informative, but are weakened by a complete lack of references. The reader who 
wishes to learn more about Positronium or the Zeroth Law of Thermodynamics is 
given no assistance here. A number of the mathematical entries, on the other hand, 
are so brief that important qualifications and clarifications are omitted. Thus, in 
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“Continuous Function” we should have 6 > 0; in “Contour Integration”’ the quan- 
tities »; should be defined; in “Analytic Function” it does not suffice for the deriva- 
tive to be single-valued at the point itself; in “Gibbs Phenomenon,” z = 2/(2n + 1) 
is not the discontinuity, but is merely near it; and in “Asymptotic Series” the ex- 
pansion may be different in different sectors in the complex plane. 

The typography is good and there appear to be relatively few errors such as 
Fronde Number in “Hydraulic Jump” and Hamiltonsion Theory in ‘Eikonal.” 

There are a few eccentricities. “Bigit’’ is advocated for what is now called “‘bit’’ 
in the binary system of numeration, and the number 7 is defined as the smallest 
positive time t at which the oscillator given by 


#=-z and 27(0)=1, «(0) =0 


again attains 7(t) = 0. 
In line with the current trend there are appropriate foreign language dictionaries 
in the appendix. The languages here are French, German, Spanish, and Russian. 
Without question, this volume will be a standard reference in many technical 
libraries. 


D.S. 


61[Z].J. F. Davison, Programming for Digital Computers, Gordon and Breach 
Science Publishers, New York, 1962, xi + 175 p., 22 cm. Price $6.00. 


The aim of this book is to provide an introduction to the general subject of writing 
programs, and it is written so as to be intelligible to the non-mathematician. It 
begins with a general discussion of the role and task of the programmer, assuming 
that he starts with the statement of a problem that needs to be programmed, and 
progresses to the point where routine computer operation has been achieved. 

The essential vehicle for discussing the techniques of programming is a theo- 
retical machine—TRIDEC—a 3-address decimal machine. With the aid of this 
machine and the limited set of orders, the author develops the basic concepts of 
programming up to the point where the reader has a feel for writing a simple 
routine using index register techniques and loops. There is then a brief discussion of 
a simple type of console to convey some notion of how the machine is controlled. 

Under the heading of more sophisticated techniques there is a look at symbolic 
programming, subroutines, and floating-point computation. 

Interpretive schemes and some aspects of automatic coding are then briefly men- 
tioned. For such a broad subject the treatment is necessarily sketchy and it at- 
tempts merely to give general impressions. 

Finally there is a discussion of differences among some different types of machines. 

Some of the operating concepts will seem odd to American programmers, par- 
ticularly the idea of using an endless loop as an equivalent to a halt. 

For its small size, the book gives a general appreciation of programming. In par- 
ticular, the details of TRIDEC coding are effectively presented. 

A. Srnkov 


National Security Agency 
Fort Meade, Maryland 
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62[Z].—D. 8. Evans, Digital Data, Their Derivation and Reduction for Analysis 
and Process Control, Interscience Publishers, Inc., New York, 1961, ii + 82 p., 
19 cm. Price $2.95. 


This amazingly small monograph presents, with typical English economy of 
words, a detailed introduction to the considerations involved in producing digital 
data from mechanical position analogs, such as shaft rotation or linear displace- 
ment. 

Chapter I, Incremental Measurements, is introductory in nature, and states the 
prime reasons and principles for automatic digitizing, namely, to conserve man- 
power and to avoid human failure in applications where both accuracy and speed 
become increasingly important. The relations among physical, graphical, and digital 
representations of quantity are discussed, as are the limitations of scaling with 
respect to ultimate accuracy and precision. The chief advantages, methods, and 
system considerations are even summarized. 

Chapters II and III, Digital Counting Devices and Direct Reading from Coded 
Scales, give brief, clear presentations of the several counting, direct reading, 
mechanical, and optical devices for analog-to-digital data conversion in the author’s 
experience. The codes employed, methods to avoid ambiguities in read-outs, and 
detailed characteristics of each device are presented. A valuable feature here is the 
listing of performance figures and system requirements for seyen specific digitizers. 

Chapter IV, Ancillary Equipment, introduces some of the additional hardware 
and techniques required in incorporating mechanical analog-to-digital converters 
into data systems. The final chapter, Chapter V, System Arrangement and Applica- 
tion, approaches the analog-to-digital data flow from an over-all system view, 
indicating some of the end uses to which decoders are often put, and which type 
of device is then selected. 

The author has succeeded in proving, that much valuable and easily applied 
information can be supplied in a very small volume, which includes 63 excellent 
figures and photographs within its 78 pages of text. 

H. M. Ernst 


Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, D. C. 


63[Z].—LesaREN Hier, Jr. & Leonarp M. Isaacson, Experimental Music 
Composition with an Electronic Computer, McGraw-Hill Book Co., Inc., New 
York, 1959, vii + 197 p., 24 cm. Price $6.00. 


This book is an exposition of the use of programming techniques, mathematics, 
and musical theory in the composition of music on the ILLIAC computer at the 
University of Illinois. It is based on the results of a set of experiments designed to 
determine whether high-speed computers could be used effectively to “compose” 
music and to analyze musical structures. No attempt was made to generate elec- 
tronic (synthetic) music, so the performance of the music composed was reserved 
for conventional musical instruments. Within the limited scope of their aims, the 
authors were successful not only in carrying out their experiments, but also in 
producing this neat and scholarly description of their work. 
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After presenting a chapter on the nature of the aesthetic problem in music, a 
brief account is given of recent technical and artistic developments in both experi- 
mentally composed music and electronically performed music. Here the authors 
mention the use of mathematical formulations to supply foundations for construct- 
ing “supposedly aesthetically significant art structures”; for example, G. D. 
Birkhoff’s theory of aesthetic measure and J. Schillinger’s theory of mathematical 
aesthetics, but none of these formulations are employed in their experiments. This 
background is followed by an analysis of the technical problems involved and a 
brief description of the programming operations of a digital computer; then a 
detailed account of each of the experiments is given. 

In developing their computer codes, the authors employed basic concepts from 
information theory and applied the Monte Carlo method. Random sequences of in- 
tegers were equated to notes in the (Western) musical scale and also to rhythmic 
patterns, dynamics, and playing techniques. These random integers were then 
screened by applying tests expressing various rules of composition (depending on 
the experiment) and accepted or rejected, depending on the rules in effect. For 
example, the first experiments began with the simplest rules for writing polyphonic 
music in C-major in cantus firmus settings, using only the 15 notes from low C to C 
above middle C. The second experiment produced four-part first-species counter- 
point in the same note range. Later experiments were based on the chromatic scale, 
and provided for the generation of rhythmic patterns, dynamic markings, and play- 
ing instructions for stringed instruments (legato, staccato, etc.). The last set of 
experiments produced what the authors call “Markoff chain music” of zero order 
and of higher orders. (The music is defined to be of zero order, if the generation 
of the interval J, between the notes n — 1 and n is independent of the choice of 
the interval J,_, ; it is of first order, if the choice of J,, is dependent on J,_, , etc.) 

The programming details are given for each experiment along with explanatory 
flow-charts, lists of composition rules employed, and tables summarizing the cal- 
culations performed. The results of the experiments, that is, the computer-con- 
structed music, were combined by the authors into a composition for string quartet, 
entitled ‘“The ILLIAC Suite,” which was transcribed by hand into conventional 
musical notation, and is given in the Appendix of the book. 

It should be emphasized that the authors do not advocate that computers should 
replace composers. They only show that computers can serve as useful aids to 
composers in experimenting with new methods and rules in the field of composition, 
in studying the classical forms of music (sonata, fugue, etc.), in analyzing the styles 
of specific composers, and in throwing some light on the aesthetic nature of music. 

Numerous references are given in the form of footnotes throughout the book; it 
would have been helpful if these references had been listed in the form of a bibliog- 
raphy at the end of the book. 

JOANNA Woop ScuHor 


Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, D. C. 
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64[Z].—Puiuir Morrison & Emity Morrison, Editors, Charles Babbage and His 
Calculating Engines, Dover Publications, Inc., New York, 1961, xxxviii + 
400 p., 21 cm. Price $2.00 (Paperbound). 


The editors have prefaced this convenient compilation of selected writings of 
Charles Babbage and of several relevant papers by his contemporaries with an 
excellent biographical sketch, followed by a brief history of punch cards and a 
selected bibliography of pertinent literature. 

The body of this entertaining book is made up of three parts, consisting, respec- 
tively, of unabridged chapters from Babbage’s autobiographical Passages from the 
Life of a Philosopher, originally published in London in 1864; essays extracted 
verbatim from Babbage’s Calculating Engines, published in London in 1889 under 
the editorship of his son Henry P. Babbage; and miscellaneous papers. These last 
include a complete listing of Babbage’s published papers, a general plan of his 
Analytical Engine, several plates relating to an eight-day clock and a hydraulic 
ram as exemplifying his notation for expressing the action of machinery, and a 
reproduction of the table of contents of Passages from the Life of a Philosopher, 
which is revelatory of certain interesting aspects of his personal life and of his multi- 
farious scientific interests. His versatility is revealed by these fields of interest, 
which include mathematics, railway engineering, cryptanalysis, and submarine 
navigation, in addition to his lifelong preoccupation with ealculating machines. 
Furthermore, his only significant completed work, namely, the book, Economy of 
Manufactures and Machinery, published in 1832, foreshadowed the field now known 
as operations research. 

The essays from Babbage’s Calculating Engines include Dr. Lardner’s detailed 
description of Babbage’s Difference Engine, originally published in 1834; a memoir 
on the Analytical Engine published by an Italian engineering officer, L. F. 
Menabrea, following a visit of Babbage to Turin in 1840, and subsequently trans- 
lated into English, with extensive annotations, by the Countess of Lovelace, only 
child of Lord Byron. The Countess, we are told, thoroughly understood and appre- 
ciated Babbage’s elaborate plans for a universal automatic digital computer, and 
is reputed to have written the most lucid contemporary accounts of that machine 
as designed by its inventor. One of her contributions, which was appended to the 
translation of Menabrea’s paper, is a detailed program for evaluating the Bernoulli 
numbers on the proposed Analytical Engine. 

In this book one reads of Babbage’s ambitious plans for an increasingly complex 
series of calculating machines, of his subsequent difficulty in securing continued 
financial support of the British Government in this research, and of his ultimate 
failure to bring his farsighted plans to fruition because of their magnitude and the 
resulting extraordinary demands on the technology of his time. 

As the editors aptly remark of Babbage in concluding their Introduction, ““The 
wide range of his practical and scientific interests and his clear commitment to the 
notion that careful analysis, mathematical procedures, and statistical calculations 
could be reliable guides in almost all facets of practical and productive life give him 
still a wonderful modernity. ... His monument, not wholly beautiful, but very 
grand, is the kind of coupled research and development that is epitomized today, as 
it was foreshadowed in his time, by the big digital computers.” 


J. W. W. 
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65[Z].—Paut Stecet, Understanding Digital Computers, John Wiley & Sons, Inc., 
New York, 1961, ix + 403 p., 23 cm. Price $8.50. 


This book is intended as a comprehensive introduction to digital computers and 
data processors for a reader with no previous knowledge of the field. The author 
does expect the reader to have a basic understanding of electronics such as that 
acquired by a “ham.” 

The material is presented in an introductory chapter and three sections: Section I 
is devoted to the basic logical elements; Section II summarizes the circuits which 
can be used as computer building blocks; Section III describes the functional parts 
of a computer, and culminates with the detailed description of a specimen com- 
puter and its use. 

The author states that the introductory chapter, which describes the digital 
computer and its functions, “provides the reader with an appreciation—a feel—for 
digital computers and a craving to learn more about them.” This reviewer is not as 
certain as the author seems to be of the reader’s reaction, for he finds the style of 
the book condescending and the material presented in a superficial fashion. 

The first example of a superficial discussion is the following: In the introduction 
addition is described in terms of counting, and multiplication is described in terms 
of repeated addition; however, in the second chapter addition tables and multiplica- 
tion tables are given for binary numbers. The relation between the two descriptions 
of addition and multiplication is never given. 

Another example has to do with complements of numbers which are discussed in 
the first section, where subtraction is done by addition of complements. However, 
the multiplication of negative numbers by means of the products of their comple- 
ments is never treated or even mentioned. Thus, the use of signed absolute values 
for the representation of numbers in a computer is never motivated. and the dis- 
cussion of multiplication in the specimen computer described in the final chapter is 
very unclear. Therein computer subtraction is carried out by use of complements 
and multiplication is done by repeated addition, but the numbers are represented 
in absolute-value form. 

The final example is furnished by the author’s discussion of instruction modifica- 
tion on pages 338 and 339. In this discussion the author states that am add-address 
instruction may be used “together with one add instruction to cause the computer 
to accumulate a long list of numbers.” The reader is never warned about the con- 
fusion that may result when the same accumulator is used for address modification 
and the accumulation of the sum, nor is he told of the price that must be paid in 
red-tape orders. Thus, he has no notion of the price that may have to be paid in 
orders and speed of computation when address modification is used. 

These examples convince the reviewer that the book under discussion does not 
answer the great need for a good introductory book on digital computers. 


a. ER. 'T. 











TABLE ERRATA 


315.—A. J. C. Cunnincuam & H. J. Woopauu, Factorization of (y” * 1), Hodgson, 
London, 1925. 
On page 6, the entry, which is given as a prime 


nm = 118, vis. 576,460,753,377,165313 
should read 5,521,693 -104,399,276,341. 
JOHN BRILLHART 


University of San Francisco 
San Francisco, California 


316.—K. M. Howe, Revised Tables of 6j-Symbols, University of Southampton, 
Mathematics Department, Research Report 59-1, March 1959. [See Review 3, 
Math. Comp., v. 14, 1960, p. 76-77.] 

A number (twenty-three) of errors have been found by L. Silverberg of Lund 
University, Sweden in these tables. The errors occurred as a result of failure to 
insert the factor 107 -107 in the printed numerator; that is, in multiplying the actual 
value of the 6j-Symbol by 107. The corrected values may be obtained by writing 
to the author. 

H. P. 


317.—D. R. Kaprexar, Cycles of Recurring Decimals, v. Il (From N = 167 to 213 
and many other numbers.) Khare Wada, Deolali, India, 1953. Published by the 
author. [See Review 1205, MT AC, v. 8, 1954, p. 148; also Table Erratum 264, 
MTAC, v. 12, 1958, p. 164-165.] 

On page 36, in the table of exponents of 10 (mod p), corresponding to P = 797, 

for C = 2, read C = 4. 

This error appears also in a table of Kraitchik [1], of which the table under dis- 
cussion is a reprint. 
ANDRZEJ MAKOwsKI 

Institute of Mathematics 


University of Warsaw 
Warsaw, Poland 


1. M. Krarroutix, Recherches sur la Théorie des Nombres, v. 1, Paris, 1924, p. 131-145. 


318.—G. W. Spencetty & R. M. Spencetzy, Smithsonian Elliptic Functions 
Tables, Smithsonian Institution, Washington, D. C., 1947. [See Review 485, 
MTAC, v. 3, 1948-1949, p. 89-92.] i 
Values of the complete elliptic integrals K and E and of Jacobi’s nome g, com- 

puted correct to 25D by the method of A. V. Hershey and the writer [1], have been 

compared with corresponding entries in the Smithsonian tables of elliptic functions. 

Errors thereby revealed in these tables can be rectified through the following correc- 

tions. 

The last (sixteenth) significant figure should be increased by a unit in the fol- 
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lowing tabular entries: 
K, for 25°; K’, for 8°; E, for 27°; E, for 29°; 
E’, for @ = 4°, 6°, 10°, 22°, 30°, 32°, 42°: 
q, for @ = 17°, 20°, 24°, 26°, 37°; 
q’, for @ = 8°, 34°, 38°, 40°. 
The last significant figure should be decreased by a unit in the following entries: 
K’', for 6 = 24°, 39°, 41°; 
E, for @ = 22°, 24°, 25°, 41°; 
q’, for @ = 11°, 28°, 29°, 31°, 37°, 39°, 44°. 
Last-figure additive corrections of greater magnitude are as follows: 
K’(3°), +5; K’(6°), —2; K’(9°), —2; 
K’(11°), —3; E’(11°), —2; ¢(3°), —19; 
q(5°), —3; q(6°), +28; g(9°), +6; 
q(10°), +2; g(11°), +5; q(15°), +2; 
q(16°), +5; 9(19°), +4; q(21°), +4; 
q(22°), +11; q(39°), +2; q(41°), +2; (44°), +2; 
q (3°), +4; (6°), —2; 9’(9°), —2; 
q'(41°), —2. 


The six errors in g corresponding to @ = 3°, 5°, 6°, 9°, 10°, and 11°, respectively, 
were originally publicized by A. Fletcher [2]. The error in g for @ = 15° appears to 
have been discovered initially by T. H. Southard and H. O. Rosay [3]. 

More serious computational errors in these tables have been discovered and cor- 
rected by G. W. Spenceley [4] and by C. R. Sexton [5]. 


A. R. DrDonato 
U. 8. Naval Weapons Laboratory 
Dahlgren, Virginia 


1. A. R. DiDonato & A. V. Hersuey, ‘‘New formulas for computing incomplete elliptic 
integrals of the first and second kind,’’ J. Assoc. Comput. Mach., v. 6, 1959, p. 515-526. 

2. A. Fietcuer, “Guide to tables of elliptic functions ” MTAC, v. 3, 1948-1040, ap 229-281. 

3. T. H. Sourarp & H. O. Rosay, MTE 281, Math. Comp. , v. 14, 1960, p. 220- 

4. G. W. Spencetny, MTE 226, MTAC, v. 7, 1953, p. 106-107. 

5. C. R. Sexton, MTE 295, Math. Comp., v. 14, 1960, p. 404. 
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CORRIGENDUM 


Dante. SHanxs, “On numbers of the form n* + 1,” Math. Comp., v. 15, 1961, 
p. 186-189. 

The entries in Table 1, on page 187, were based upon an erroneous count of 
primes of the form n* + 1. Two primes, corresponding to n = 786 and n = 788, 
respectively, were omitted; they are listed in a recent paper [1] of Gloden. The 
necessary changes in Table 1 are as follows: Q,(N) = 96, 100, and 111 for N = 800, 
900, and 1000, respectively. The corresponding values of Q,/Q, , should then read 
0.974, 0.922, and 0.938. The agreement of the proposed asymptotic formula for 
Q,(N) is, therefore, improved. On the second line of section 2, “The 109th prime” 
should be changed to “The 111th prime.” 

D. 8. 


1. A. GLopEn, “Additions to Cunningham’s factor table of N‘ + 1,”” Math. Comp., v. 16, 
1962, p. 239-241. 
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